1-Notas aclaratorias y breve resumen

BiciMAD es el sistema de bicicletas electricas de Madrid. Es uno de los medios mas novedosos y limpios de la capital y podria suponer una autentica revolucion en el futuro de la movilidad urbana. A 31 de Enero de 2020, cuenta con mas de 62000 abonados anuales y lo han usado hasta 148000 usuarios ocasionales desde el inicio del sistema.

Este archivo es un breve analisis de las cuestiones mas importantes que rodean a BiciMAD. En el mismo, se han omitido caracteres especiales del castellano como tildes, apostrofes etc, debido al encoding UTF-8. Disculpen las molestias.

Los datos se han obtenido del portal de datos abiertos de Madrid, OPENDATA, subidos por la EMT de Madrid y se pueden encontrar en la siguiente URL: https://opendata.emtmadrid.es/Datos-estaticos/Datos-generales-(1)

2-Creacion del archivo a analizar

Lectura de nuestros archivos mensuales CSV, una vez convertidos. Leeremos los archivos mensuales uno a uno, agregando un separador “;” para que se lea correctamente, en caso de necesitarlo.

#Enero<-read.csv("201801_BICIMAD_GEO.csv", sep=";")
#Febrero<-read.csv("201802_BICIMAD_GEO.csv", sep=";")
#Marzo<-read.csv("201803_BICIMAD_GEO.csv", sep=";")
#Abril<-read.csv("201804_BICIMAD_GEO.csv", sep=";")
#Mayo<-read.csv("201805_BICIMAD_GEO.csv", sep=";")
#Junio<-read.csv("201806_BICIMAD_GEO.csv", sep=";")
#Julio<-read.csv("201807_BICIMAD_GEO.csv", sep=";")
#Agosto<-read.csv("201808_BICIMAD_GEO.csv", sep=";")
#Septiembre<-read.csv("201809_BICIMAD_GEO.csv", sep=";")
#Octubre<-read.csv("201810_BICIMAD_GEO.csv", sep=";")
#Noviembre<-read.csv("201811_BICIMAD_GEO.csv", sep=";")
#Diciembre<-read.csv("201812_BICIMAD_GEO.csv",sep=";")

2.1-Juntaremos los archivos

Con la funcion plyr::rbind.fill, que juntara las filas de nuestros archivos mensuales en una gran base de datos

*Nota:Tambien podemos realizar el mismo proceso con el paquete dplyr y la funcion bind_rows o rbind utilizando base R

#library(plyr)
#library(utils)

#Juntamos los archivos en un solo "dataset" para trabajar mas facilmente.

#BiciMAD<-rbind.fill(Enero,Febrero,Marzo,Abril,Mayo,Junio,Julio,Agosto,Septiembre,Octubre,Noviembre,Diciembre)

#Para liberar memoria RAM, quitaremos los archivos mensuales del Global Environment

#rm(Enero)
#rm(Febrero)
#rm(Marzo)
#rm(Abril)
#rm(Mayo)
#rm(Junio)
#rm(Julio)
#rm(Agosto)
#rm(Septiembre)
#rm(Octubre)
#rm(Noviembre)
#rm(Diciembre)

##Una vez realizado esto, podemos guardar el archivo resultante "BiciMAD" como CSV para cargarlo directamente cuando procedamos al analisis 
#utils::write.csv(BiciMAD,"BiciMAD.csv")

#quitaremos el paquete plyr puesto que puede causar problemas con el paquete dplyr que posteriormente utilizaremos.
#detach("package:plyr")

3-Primeros comandos

#Ahora, abriremos el archivo BiciMAD que hemos creado
BiciMAD<-read.csv("BiciMAD.csv")

#Comenzaremos mirando la estructura de nuestros datos
str(BiciMAD)
## 'data.frame':    3678086 obs. of  19 variables:
##  $ X               : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ X_id            : Factor w/ 3678086 levels "5a505acd2f384319304ed614",..: 1 2 3 4 6 10 5 7 8 9 ...
##  $ ageRange        : int  5 0 0 3 4 5 0 4 4 5 ...
##  $ idplug_base     : int  1 4 19 1 5 5 21 20 1 12 ...
##  $ idplug_station  : int  7 117 117 110 58 74 18 17 63 134 ...
##  $ address_end     : Factor w/ 173 levels "#N/A","Avenida de Alberto Alcocer nº 22",..: 59 23 23 109 30 7 144 98 49 128 ...
##  $ long_end        : Factor w/ 173 levels "-3.6689089","-3.6691523",..: 97 167 167 47 106 23 127 143 130 136 ...
##  $ lat_end         : Factor w/ 173 levels "#N/A","40.3919385",..: 95 139 139 112 111 49 90 88 104 13 ...
##  $ idunplug_base   : int  14 21 19 2 3 18 1 7 16 23 ...
##  $ idunplug_station: int  6 24 24 82 169 92 34 128 3 135 ...
##  $ address_start   : Factor w/ 173 levels "#N/A","Avenida de Alberto Alcocer nº 22",..: 53 82 82 35 72 29 91 90 81 127 ...
##  $ long_start      : Factor w/ 173 levels "-3.6689089","-3.6691523",..: 114 111 111 14 103 36 154 79 132 91 ...
##  $ lat_start       : Factor w/ 173 levels "#N/A","40.3919385",..: 123 69 69 7 132 119 67 10 74 2 ...
##  $ travel_time     : int  284 666 662 708 171 427 166 722 254 305 ...
##  $ unplug_hourDate : Factor w/ 358 levels "01/01/2018","01/01/2019",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ unplug_hourTime : Factor w/ 24 levels "0:00:00","1:00:00",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ user_day_code   : Factor w/ 2049401 levels "00001f08f2d92c952bbb1a3f69926df6cd2fc414d6fa6c999b15f3dd406e1470",..: 1409169 812935 653040 87086 639772 514909 1898701 19003 208312 1465617 ...
##  $ user_type       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ zip_code        : int  28010 NA NA 28005 28004 28006 NA 28005 28004 28012 ...
#Como vemos, algunas variables estan codificadas con clases distintas a lo que deberian ser. Procedemos a corregir estas variables

#Factores
BiciMAD$ageRange<-as.factor(BiciMAD$ageRange)
BiciMAD$idplug_base<-as.factor(BiciMAD$idplug_base)
BiciMAD$idplug_station<-as.factor(BiciMAD$idplug_station)
BiciMAD$idunplug_base<-as.factor(BiciMAD$idunplug_base)
BiciMAD$idunplug_station<-as.factor(BiciMAD$idunplug_station)

#Codificaremos las variables de latitud y longitud para ver si sus niveles corresponden con address_start y address_end
BiciMAD$long_end<-as.factor(BiciMAD$long_end)
BiciMAD$long_start<-as.factor(BiciMAD$long_start) 
BiciMAD$lat_end<-as.factor(BiciMAD$lat_end)
BiciMAD$lat_start<-as.factor(BiciMAD$lat_start)

BiciMAD$user_type<-as.factor(BiciMAD$user_type)
BiciMAD$zip_code<-as.factor(BiciMAD$zip_code)

#Fecha de salida de la bici, especificando formato dia/mes/ano
BiciMAD$unplug_hourDate<-as.Date(BiciMAD$unplug_hourDate, format="%d/%m/%Y")

#Quitaremos variables que NO nos dan informacion util (X_id y user_day_code)
BiciMAD$X_id<-NULL
BiciMAD$user_day_code<-NULL

#Variables numericas
BiciMAD$travel_time<-as.numeric(BiciMAD$travel_time)

A priori, parece que en este dataset son todo factores o variables categoricas. Al ser un dataset que nos indica los puntos de A a B, muchas variables cuentan con los mismos casos(Ej. si un individuo u otro dejan la bici en Plaza del sol, las variables de address, longitud y latitud seran las mismas tanto para uno como para el otro.)

La unica variable que no se repite y que es medida (independiente para cada usuario) es travel_time que mide el tiempo en segundos de A a B. Tambien he ajustado la variable de fecha como fecha y no como factor.

3.1-Correccion de errores/NAs variable a variable

Procedemos a explorar estadistica simple y si hay NA’s en nuestra base de datos. Realizamos un “str” y un “summary” para comprobar cambios.

#Funcion de estrutura
str(BiciMAD)
## 'data.frame':    3678086 obs. of  17 variables:
##  $ X               : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ ageRange        : Factor w/ 7 levels "0","1","2","3",..: 6 1 1 4 5 6 1 5 5 6 ...
##  $ idplug_base     : Factor w/ 30 levels "1","2","3","4",..: 1 4 19 1 5 5 21 20 1 12 ...
##  $ idplug_station  : Factor w/ 173 levels "1","2","3","4",..: 7 114 114 107 57 71 18 17 62 131 ...
##  $ address_end     : Factor w/ 173 levels "#N/A","Avenida de Alberto Alcocer nº 22",..: 59 23 23 109 30 7 144 98 49 128 ...
##  $ long_end        : Factor w/ 173 levels "-3.6689089","-3.6691523",..: 97 167 167 47 106 23 127 143 130 136 ...
##  $ lat_end         : Factor w/ 173 levels "#N/A","40.3919385",..: 95 139 139 112 111 49 90 88 104 13 ...
##  $ idunplug_base   : Factor w/ 30 levels "1","2","3","4",..: 14 21 19 2 3 18 1 7 16 23 ...
##  $ idunplug_station: Factor w/ 173 levels "1","2","3","4",..: 6 23 23 79 166 89 33 125 3 132 ...
##  $ address_start   : Factor w/ 173 levels "#N/A","Avenida de Alberto Alcocer nº 22",..: 53 82 82 35 72 29 91 90 81 127 ...
##  $ long_start      : Factor w/ 173 levels "-3.6689089","-3.6691523",..: 114 111 111 14 103 36 154 79 132 91 ...
##  $ lat_start       : Factor w/ 173 levels "#N/A","40.3919385",..: 123 69 69 7 132 119 67 10 74 2 ...
##  $ travel_time     : num  284 666 662 708 171 427 166 722 254 305 ...
##  $ unplug_hourDate : Date, format: "2018-01-01" "2018-01-01" ...
##  $ unplug_hourTime : Factor w/ 24 levels "0:00:00","1:00:00",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ user_type       : Factor w/ 4 levels "0","1","2","3": 2 2 2 2 2 2 2 2 2 2 ...
##  $ zip_code        : Factor w/ 768 levels "0","2","28","112",..: 267 NA NA 262 261 263 NA 262 261 269 ...
#Funcion de summary
summary(BiciMAD)
##        X           ageRange     idplug_base      idplug_station   
##  Min.   :      1   0:1355749   1      : 248458   135    :  53898  
##  1st Qu.: 919522   1:  46335   2      : 222518   175    :  51297  
##  Median :1839044   2:  41848   3      : 208513   43     :  51242  
##  Mean   :1839044   3: 236331   4      : 186692   129    :  50132  
##  3rd Qu.:2758565   4:1217575   24     : 182486   163    :  49562  
##  Max.   :3678086   5: 763320   5      : 171424   57     :  48394  
##                    6:  16928   (Other):2457995   (Other):3373561  
##                                    address_end            long_end      
##  Paseo de la Chopera nº 14               :  53898   -3.6971829:  53898  
##  Calle Jaime el Conquistador nº30        :  51297   -3.6981528:  51297  
##  Plaza de la Cebada nº 16                :  51242   -3.7088337:  51242  
##  Paseo de Santa María de la Cabeza nº 58 :  50132   -3.6987665:  50132  
##  Calle del General Álvarez de Castro nº 2:  49562   -3.7015686:  49562  
##  Calle Valencia nº 1                     :  48394   -3.7008803:  48394  
##  (Other)                                 :3373561   (Other)   :3373561  
##        lat_end        idunplug_base     idunplug_station 
##  40.3919385:  53898   1      : 247966   135    :  52550  
##  40.3961981:  51297   2      : 221545   175    :  50888  
##  40.4112744:  51242   3      : 207771   43     :  50659  
##  40.4017926:  50132   4      : 186420   129    :  49472  
##  40.4344731:  49562   24     : 180642   163    :  49366  
##  40.4089282:  48394   5      : 170990   57     :  48673  
##  (Other)   :3373561   (Other):2462752   (Other):3376478  
##                                   address_start          long_start     
##  Paseo de la Chopera nº 14               :  52550   -3.6971829:  52550  
##  Calle Jaime el Conquistador nº30        :  50888   -3.6981528:  50888  
##  Plaza de la Cebada nº 16                :  50659   -3.7088337:  50659  
##  Paseo de Santa María de la Cabeza nº 58 :  49472   -3.6987665:  49472  
##  Calle del General Álvarez de Castro nº 2:  49366   -3.7015686:  49366  
##  Calle Valencia nº 1                     :  48673   -3.7008803:  48673  
##  (Other)                                 :3376478   (Other)   :3376478  
##       lat_start        travel_time      unplug_hourDate      unplug_hourTime   
##  40.3919385:  52550   Min.   :      0   Min.   :2018-01-01   19:00:00: 270285  
##  40.3961981:  50888   1st Qu.:    453   1st Qu.:2018-04-22   18:00:00: 260738  
##  40.4112744:  50659   Median :    689   Median :2018-07-17   20:00:00: 245184  
##  40.4017926:  49472   Mean   :   1131   Mean   :2018-07-11   15:00:00: 235315  
##  40.4344731:  49366   3rd Qu.:   1061   3rd Qu.:2018-10-03   8:00:00 : 234449  
##  40.4089282:  48673   Max.   :5730890   Max.   :2019-01-01   14:00:00: 231202  
##  (Other)   :3376478                                          (Other) :2200913  
##  user_type      zip_code      
##  0:   2097   28005  : 215143  
##  1:3327888   28012  : 191429  
##  2:  59990   28007  : 175903  
##  3: 288111   28004  : 162836  
##              28009  : 134675  
##              (Other):1200192  
##              NA's   :1597908
#Comprobar NAs en variables de perfil del usuario
table(is.na(BiciMAD$ageRange))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$user_type))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$zip_code))
## 
##   FALSE    TRUE 
## 2080178 1597908
#Comprobar NAs en variables ORIGEN DEL RECORRIDO
table(is.na(BiciMAD$idunplug_base))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$idunplug_station))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$address_start))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$long_start))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$lat_start))
## 
##   FALSE 
## 3678086
#Comprobar NAs en variables DESTINO DEL RECORRIDO
table(is.na(BiciMAD$idplug_base))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$idplug_station))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$address_end))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$long_end))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$lat_end))
## 
##   FALSE 
## 3678086
#Comprobar NAs en otras variables de tiempo
table(is.na(BiciMAD$travel_time))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$unplug_hourDate))
## 
##   FALSE 
## 3678086
table(is.na(BiciMAD$unplug_hourTime))
## 
##   FALSE 
## 3678086

Como podemos observar, la unica variable con NAs es “zip_code”. Recordemos que anteriormente hicimos un VLOOKUP/INDEXMATCH en EXCEL para crear las variables de “adress”, “long” y “lat”. Si no encontraba el registro de la estacion, Excel devolvia #N/A! (escrito y no un NA como variable nula) por lo que debemos buscar estos valores en nuestra base de datos y eliminarlos.

En nuestro primer str() hemos podido ver que address tiene 173 niveles como factor, lo mismo que en long y lat. Sin embargo, cuando realice el VLOOKUP/INDEXMATCH en Excel, solo tenia 172 estaciones oficialmente (tanto en el archivo de bases que relacionaba estaciones con IDs como los archivos de las bases de datos con la variable idplug y unplug con las que estoy trabajando). Esto significa por lo tanto, que hay UNA estacion que sobra.

#Este comando creara otro dataframe donde podemos ver que efectivamente, hay otra estacion oculta en los datos. La numero 2008.
NA_escrito<-subset(BiciMAD, address_start=="#N/A") 

summary(NA_escrito)
##        X           ageRange  idplug_base idplug_station
##  Min.   : 131701   0:1      1      :8    2008   :17    
##  1st Qu.:2620711   1:9      2      :8    1      : 0    
##  Median :2995852   2:0      3      :1    2      : 0    
##  Mean   :2582820   3:0      4      :0    3      : 0    
##  3rd Qu.:3021959   4:7      5      :0    4      : 0    
##  Max.   :3563294   5:0      6      :0    5      : 0    
##                    6:0      (Other):0    (Other): 0    
##                            address_end       long_end        lat_end  
##  #N/A                            :17   #N/A      :17   #N/A      :17  
##  Avenida de Alberto Alcocer nº 22: 0   -3.6689089: 0   40.3919385: 0  
##  Avenida de Alfonso XII nº 54    : 0   -3.6691523: 0   40.3961981: 0  
##  Avenida de Menéndez Pelayo nº 11: 0   -3.6691838: 0   40.3972616: 0  
##  Avenida de Menéndez Pelayo nº 3 : 0   -3.6697526: 0   40.3974119: 0  
##  Avenida de Menéndez Pelayo nº 38: 0   -3.670422 : 0   40.400781 : 0  
##  (Other)                         : 0   (Other)   : 0   (Other)   : 0  
##  idunplug_base idunplug_station                          address_start
##  1      :9     2008   :17       #N/A                            :17   
##  2      :7     1      : 0       Avenida de Alberto Alcocer nº 22: 0   
##  3      :1     2      : 0       Avenida de Alfonso XII nº 54    : 0   
##  4      :0     3      : 0       Avenida de Menéndez Pelayo nº 11: 0   
##  5      :0     4      : 0       Avenida de Menéndez Pelayo nº 3 : 0   
##  6      :0     5      : 0       Avenida de Menéndez Pelayo nº 38: 0   
##  (Other):0     (Other): 0       (Other)                         : 0   
##       long_start      lat_start   travel_time      unplug_hourDate     
##  #N/A      :17   #N/A      :17   Min.   :      9   Min.   :2018-01-18  
##  -3.6689089: 0   40.3919385: 0   1st Qu.:    359   1st Qu.:2018-09-24  
##  -3.6691523: 0   40.3961981: 0   Median :    805   Median :2018-10-22  
##  -3.6691838: 0   40.3972616: 0   Mean   : 196432   Mean   :2018-09-16  
##  -3.6697526: 0   40.3974119: 0   3rd Qu.:   3879   3rd Qu.:2018-10-24  
##  -3.670422 : 0   40.400781 : 0   Max.   :3182807   Max.   :2018-12-17  
##  (Other)   : 0   (Other)   : 0                                         
##  unplug_hourTime user_type    zip_code 
##  16:00:00:5      0: 0      0      : 0  
##  9:00:00 :5      1:10      2      : 0  
##  11:00:00:3      2: 0      28     : 0  
##  10:00:00:1      3: 7      112    : 0  
##  14:00:00:1                135    : 0  
##  15:00:00:1                (Other): 0  
##  (Other) :1                NA's   :17
#Con el comando summary, una vez mas, podemos ver que estos registros se dan en el mes de Octubre, con mas casos sobre las 9:00 y 16:00, donde las bases mas utilizadas para desengancharlas son la primera, segunda y tercera y el tipo de usuario suele ser 1(suscriptor) o 3(empleado). Es posible que esto sea un almacen o bien un simple fallo, por lo que procedere a borrar todos los registros de la estacion 2008

rm(NA_escrito)

BiciMAD_sin2008<-subset(BiciMAD, address_start!="#N/A") 
str(BiciMAD_sin2008)
## 'data.frame':    3678069 obs. of  17 variables:
##  $ X               : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ ageRange        : Factor w/ 7 levels "0","1","2","3",..: 6 1 1 4 5 6 1 5 5 6 ...
##  $ idplug_base     : Factor w/ 30 levels "1","2","3","4",..: 1 4 19 1 5 5 21 20 1 12 ...
##  $ idplug_station  : Factor w/ 173 levels "1","2","3","4",..: 7 114 114 107 57 71 18 17 62 131 ...
##  $ address_end     : Factor w/ 173 levels "#N/A","Avenida de Alberto Alcocer nº 22",..: 59 23 23 109 30 7 144 98 49 128 ...
##  $ long_end        : Factor w/ 173 levels "-3.6689089","-3.6691523",..: 97 167 167 47 106 23 127 143 130 136 ...
##  $ lat_end         : Factor w/ 173 levels "#N/A","40.3919385",..: 95 139 139 112 111 49 90 88 104 13 ...
##  $ idunplug_base   : Factor w/ 30 levels "1","2","3","4",..: 14 21 19 2 3 18 1 7 16 23 ...
##  $ idunplug_station: Factor w/ 173 levels "1","2","3","4",..: 6 23 23 79 166 89 33 125 3 132 ...
##  $ address_start   : Factor w/ 173 levels "#N/A","Avenida de Alberto Alcocer nº 22",..: 53 82 82 35 72 29 91 90 81 127 ...
##  $ long_start      : Factor w/ 173 levels "-3.6689089","-3.6691523",..: 114 111 111 14 103 36 154 79 132 91 ...
##  $ lat_start       : Factor w/ 173 levels "#N/A","40.3919385",..: 123 69 69 7 132 119 67 10 74 2 ...
##  $ travel_time     : num  284 666 662 708 171 427 166 722 254 305 ...
##  $ unplug_hourDate : Date, format: "2018-01-01" "2018-01-01" ...
##  $ unplug_hourTime : Factor w/ 24 levels "0:00:00","1:00:00",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ user_type       : Factor w/ 4 levels "0","1","2","3": 2 2 2 2 2 2 2 2 2 2 ...
##  $ zip_code        : Factor w/ 768 levels "0","2","28","112",..: 267 NA NA 262 261 263 NA 262 261 269 ...

Antes de comenzar a trabajar con los datos realizaremos otro ajuste importante. Hay veces en las que un usuario por desconocimiento del sistema de BiciMAD puede desanclar y volver a anclar una bicicleta en el mismo sitio Y DONDE LA VARIABLE travel_time SUELE SER MUY PEQUENA, esto es un falso positivo. Hay otras veces donde desanclar y anclar de nuevo la bicicleta con un travel_time elevado puede significar que el usuario ha ido a hacer algo y ha vuelto, probablemente porque viva cerca de la direccion de origen/destino.

Por ello, voy a hacer un “subset” donde las estaciones “unplug” y “plug” (origen y destino) no sean la misma O BIEN, donde la variable de tiempo “travel_time” supere los 90 segundos, un minuto y medio, porque puede ser factible realizar el recorrido si las estaciones estan muy cercanas en este tiempo. De la misma manera, tambien es factible que alguien haga un recorrido largo donde origen y destino sean iguales. Una vez realizado esto, podemos observar que hay otros viajes donde la estacion de origen y destino no es la misma y el travel_time es inferior a los 45 segundos, esto puede deberse a fallos en el sistema de BiciMAD, por lo que no borrare dichos registros.

#Eliminacion de falsos positivos (misma estacion plug y unplug o trayectos con menos de 90 segundos de recorrido)
BiciMAD_sin2008_travel_timeadj<-subset(BiciMAD_sin2008, idunplug_station!=idplug_station | travel_time>=90)

rm(BiciMAD)
rm(BiciMAD_sin2008)

4-Analisis exploratorio de BiciMAD

4.1-Ajustes preliminares

#Cargaremos la mayoria de librerias a utilizar
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(scales)
library(geosphere)
library(shiny)
## 
## Attaching package: 'shiny'
## The following object is masked from 'package:geosphere':
## 
##     span
library(Rcpp)
library(leaflet)
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(tidyverse)
## -- Attaching packages --------------------------------------------------------------------------------- tidyverse 1.3.0 --
## v tibble  3.0.1     v purrr   0.3.4
## v tidyr   1.1.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## -- Conflicts ------------------------------------------------------------------------------------ tidyverse_conflicts() --
## x readr::col_factor() masks scales::col_factor()
## x purrr::discard()    masks scales::discard()
## x dplyr::filter()     masks stats::filter()
## x dplyr::lag()        masks stats::lag()
#Reasignamos el dataset ajustado a un nombre mas corto, BiciMAD
BiciMAD<-BiciMAD_sin2008_travel_timeadj

#Borraremos el antiguo dataset para liberar memoria RAM
rm(BiciMAD_sin2008_travel_timeadj)

4.2-Analisis exploratorio de la variable medible, tiempo de recorrido o travel_time

min(BiciMAD$travel_time) 
## [1] 0
max(BiciMAD$travel_time)
## [1] 5730890
#Como podemos observar el minimo registro en la variable es 0. Esto puede deberse a fallos de medicion en la bicicleta, de la misma manera el maximo es de 5730890 segundos, decidimos investigarlo.

head(arrange(BiciMAD,desc(travel_time)), n = 10)
##          X ageRange idplug_base idplug_station                   address_end
## 1  2752110        5           5              4 Plaza del Conde Suchil nº 2-4
## 2  2733274        0           9            134    Paseo de la Esperanza nº 2
## 3  1653343        5           8            151 Paseo de la Castellana nº 106
## 4   122239        1          17            108    Calle Ortega y Gasset nº 6
## 5  2827823        5          13            108    Calle Ortega y Gasset nº 6
## 6  3482291        4          22            154  Paseo de la Castellana nº 67
## 7  3482290        5          23            154  Paseo de la Castellana nº 67
## 8  3482289        4          24            154  Paseo de la Castellana nº 67
## 9  2785352        4          15            111       Calle Serrano nº 34 - A
## 10 2298733        0          15            139 Calle del General Yagüe nº 57
##      long_end    lat_end idunplug_base idunplug_station
## 1  -3.7069171 40.4302937            14              131
## 2  -3.7064516 40.4035988            17               42
## 3   -3.690861 40.4453307             2               67
## 4  -3.6865654 40.4303057            15              113
## 5  -3.6865654 40.4303057            24              108
## 6  -3.6917932 40.4457414            22              159
## 7  -3.6917932 40.4457414            18              159
## 8  -3.6917932 40.4457414             6              151
## 9  -3.6877227 40.4251002             7               91
## 10 -3.7009675 40.4572824            18              174
##                        address_start long_start  lat_start travel_time
## 1         Calle Guzmán el Bueno nº 2 -3.7133412 40.4306458     5730890
## 2      Carrera de San Francisco nº 1 -3.7120734 40.4110406     2833493
## 3                  Calle Ibiza nº 62 -3.6708959 40.4179237     2755282
## 4                 Calle Serrano nº 6 -3.6884369 40.4215246     1651625
## 5         Calle Ortega y Gasset nº 6 -3.6865654 40.4303057     1625256
## 6  Calle José Gutiérrez Abascal nº 2 -3.6907784 40.4396792     1603461
## 7  Calle José Gutiérrez Abascal nº 2 -3.6907784 40.4396792     1603440
## 8      Paseo de la Castellana nº 106  -3.690861 40.4453307     1603113
## 9               Paseo del Prado nº 1 -3.6926217 40.4186148     1582306
## 10              Paseo Yeserías nº 45 -3.7084067 40.3974119     1346583
##    unplug_hourDate unplug_hourTime user_type zip_code
## 1       2018-10-03         8:00:00         3    28003
## 2       2018-10-01        21:00:00         3     <NA>
## 3       2018-07-03        12:00:00         3     <NA>
## 4       2018-01-17        12:00:00         3     <NA>
## 5       2018-10-08        12:00:00         3     <NA>
## 6       2018-12-09         1:00:00         3     <NA>
## 7       2018-12-09         1:00:00         3     <NA>
## 8       2018-12-09         1:00:00         3     <NA>
## 9       2018-10-05         7:00:00         3     <NA>
## 10      2018-09-02        17:00:00         3     <NA>
#Como podemos observar los registros mas altos son los de empleados o user_type = 3, probablemente por labores de mantenimiento y recolocacion de las bicicletas. De la misma manera, podemos ver donde las recolocan, con la Calle Gutierrez Abascal numero 2 o Guzman el Bueno numero 2 como origen y Paseo de la Castellana numero 67 y Calle Ortega y Gasset numero 6 como destino. A la hora de calcular la asimetria de las estaciones, NO quitaremos ningun registro de los empleados ya que queremos ver la asimetria real despues de la labor de los mismos.

Tiempo de recorrido o travel_time y su relacion con los grupos de edad o ageRange

#Grafico rapido para ver el recuento de usuarios segun su grupo de edad. Como vemos, no todos los grupos de edad son representativos, ya que el grupo 1,2 y 6 contienen muy pocos registros. A su vez podemos ver que los registros "0" o "no asignado" son muy numerosos.

ggplot(BiciMAD, aes(ageRange))+geom_bar()+ scale_y_continuous(labels = comma)+xlab("Grupo de edad")+ylab("Numero de usuarios")

#Calculo de medias de la variable tiempo de recorrido, respecto al grupo de edad. Nos centramos en los grupos 3,4 y 5 que son los grupos con mas registros y donde podemos capturar suficiente informacion para sacar conclusiones. Por lo que vemos, los usuarios del grupo 3 realizan recorridos mas cortos que los del grupo de edad 4 y estos menos que los del grupo de edad 5, antes de sacar estas conclusiones, exploraremos en mas detalle.

tapply(BiciMAD$travel_time,BiciMAD$ageRange, mean, na.rm=T)
##         0         1         2         3         4         5         6 
## 1204.0700 3410.1417  887.1594  961.5068 1049.3673 1186.3557 1143.4194
#Ante tal diferencia en las medias computadas, decido hacer un subset conteniendo unicamente los usuarios con un pase anual, a sabiendas de que el user_type = 3 o empleado y user_type=2 u ocasional podrian distorsionar la media debido a los grandes valores de tiempo de recorrido en la variable: travel_time. 

#Decido hacer la misma media, donde anadiremos un filtro donde la duracion aproximada debe ser superior a los 30 segundos para que no afecte la media general, ya que hay multiples casos donde travel_time no ha sido medida bien. Vemos que el resultado es mas parejo, con medias mucho mas igualadas y sin una gran diferencia entre grupo de edades. 

#De todos modos, podemos decir que los usuarios del grupo 5, son uno de los grupos que mas tiempo utilizan la bicicleta por recorrido efectuado. Los usuarios del grupo de edad 2 y 6 consiguen una media elevada aunque no suponen un gran porcentaje de los registros en nuestros datos. Respecto al gran valor en el grupo de edad 6, teniendo mas de 66 anos con pase anual (user_type = 1) tener una media de tiempo de recorrido tan alta podria significar que la velocidad a la que se efectua el recorrido es mucho menor. 


Solo_paseanual<-subset(BiciMAD,user_type=="1" & travel_time>=30)

tapply(Solo_paseanual$travel_time,Solo_paseanual$ageRange, mean, na.rm=T)
##         0         1         2         3         4         5         6 
## 1005.7479  831.8858  886.6284  857.2191  809.3744  882.3469 1142.8863
#Grafico de la variable tiempo de recorrido, respecto al grupo de edad

ggplot(BiciMAD, aes(ageRange, travel_time)) + geom_violin()+ylim(0,7500)+xlab("Grupo de edad")+ylab("Duracion del recorrido (en segundos)")
## Warning: Removed 30609 rows containing non-finite values (stat_ydensity).

#Una vez mas, podemos ver que no hay diferencias notables, hemos acotado travel_time hasta los 7500 segundos, debido a los valores extremos/outliers que presenta la variable. En el grupo 1 podemos ver que no hay una gran cantidad de personas en la media, con datos muy dispersos y con valores altos y bajos.

#Limpiamos
rm(Solo_paseanual)

Tiempo de recorrido o travel_time y su relacion con Tipo de usuario o user_type

#Calculo de medias de la variable tiempo de recorrido, respecto al tipo de usuario

tapply(BiciMAD$travel_time, BiciMAD$user_type, mean, na.rm=T)
##         0         1         2         3 
## 1956.0787  903.0235 2566.2672 4306.7681
#No tendremos en cuenta los user_type = 0, pero podemos obtener conclusiones sobre los demas tipos de usuarios. El user_type = 1 o subscritor de pase anual tiene una media muy inferior al user_type = 2 o usuario ocasional/turista. Esto es logico ya que los subscriptores de pase anual cogeran mas la bici para trayectos mas cortos y de menor duracion mientras que los usuarios ocasionales cogeran la bici para trayectos mas largos y de mayor duracion. El user_type = 3 o empleado es el que mas tiempo de recorrido tiene al tener que realizar la recolocacion de bicicletas.


#Grafico de la variable tiempo de recorrido, respecto al tipo de usuario
ggplot(BiciMAD, aes(user_type, travel_time)) + geom_boxplot()+ylim(0,7500)+xlab("Tipo de usuario")+ylab("Duracion del recorrido (en segundos)")
## Warning: Removed 30609 rows containing non-finite values (stat_boxplot).

#En este boxplot, se puede observar el distinto comportamiento de los usuarios respecto a la duracion del recorrido. El tipo 1 tiene una media muy inferior a la del tipo 2 y 3. El tipo 1 tiene una simetria muy pareja con los cuartiles 1 y 3, reflejando un comportamiento mas equitativo respecto a la duracion del recorrido, mientras que la distribucion del tipo 2 es mas dispersa, debido al diferente comportamiento de los usuarios ocasionales.

4.2-Relacion entre Grupo de edad o ageRange y tipo de usuario o user_type

#Relacionamos ahora el tipo de usuario y su rango de edad. Lo pondremos en porcentaje ya que en terminos absolutos, hay muchos mas valores de 1 (subscriptor) que de 0 (sin identificar), 2 (usuario ocasional) y 3 (empleados), logicamente.

ggplot(BiciMAD, aes(x=user_type))+geom_bar(aes(fill=ageRange), position = "fill")+xlab("Tipo de usuario")+ylab("Porcentaje(%) de usuarios")

#0:No se ha podido determinar el rango de edad del usuario
#1: El usuario tiene entre 0 y 16 anos
#2: El usuario tiene entre 17 y 18 anos
#3: El usuario tiene entre 19 y 26 anos
#4: El usuario tiene entre 27 y 40 anos
#5: El usuario tiene entre 41 y 65 anos
#6: El usuario tiene 66 anos o mas 

#Como podemos ver, gran parte de los usuarios 2 u ocasionales, no se identifican con su edad al utilizar el servicio. De la misma manera, a gran parte de los usuarios 1 o subscriptores no se les puede identificar. De los identificados como usuarios 1, mas del 50% de usuarios los componen el grupo 4 y 5. El grupo 3 tambien compone alrededor de un 10% aproximadamente.

#Tambien se debe tener en cuenta que el ageRange o grupo de edad 2 lo componen usuarios de 17 y 18 anos de edad, un intervalo que agrupa unicamente 2 anos y que es muy pequeno si se compara con el ageRange o grupo de edad 4 o 5, que agrupan 13 y 24 anos de intervalo.

4.3-Creacion de mas variables

#Creacion de la variable Distancia en linea recta y Velocidad aproximada
BiciMAD$long_start<-as.numeric(as.character(BiciMAD$long_start))
BiciMAD$lat_start<-as.numeric(as.character(BiciMAD$lat_start))
BiciMAD$long_end<-as.numeric(as.character(BiciMAD$long_end))
BiciMAD$lat_end<-as.numeric(as.character(BiciMAD$lat_end))


#Calculamos la distancia en linea recta, esto tiene una serie de consecuencias
BiciMAD$SL_distancia <- distHaversine(BiciMAD[,c("long_start","lat_start")],BiciMAD[,c("long_end","lat_end")])
#Calculamos la variable distancia en linea recta entre dos estaciones. Logicamente la distancia real sera algo mayor. Tambien debemos tener en cuenta que en los casos donde las estaciones donde se ancla y se desancla la bicicleta sea la misma, la funcion dara un 0/cero.

BiciMAD$Velocidad<- BiciMAD$SL_distancia/BiciMAD$travel_time #en metros/segundo
#Podremos crear la variable velocidad, ya que tenemos la distancia en linea recta y el tiempo del viaje. Esta Velocidad tiene fallos por lo explicado anteriormente. La variable travel_time no esta bien contabilizada tampoco en muchos casos.

#Distancia entre estaciones respecto a grupo de edad
ggplot(BiciMAD, aes(x=ageRange, y=SL_distancia))+geom_violin()+xlab("Grupo de edad")+ylab("Distancia (en metros)")+scale_y_continuous(labels=comma)

#Distancia promedio respecto a tipo de usuario
ggplot(BiciMAD, aes(x=user_type, y=SL_distancia))+geom_violin()+xlab("Tipo de usuario")+ylab("Distancia (en metros)")+scale_y_continuous(labels=comma)

#Velocidad promedio respecto a grupo de edad
ggplot(BiciMAD, aes(x=ageRange, y=Velocidad))+geom_violin()+xlab("Grupo de edad")+ylab("Velocidad (en m/s)")+scale_y_continuous(labels=comma)+ylim(0,15)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 583 rows containing non-finite values (stat_ydensity).

#Velocidad promedio respecto a tipo de usuario
ggplot(BiciMAD, aes(x=user_type, y=Velocidad))+geom_violin()+xlab("Tipo de usuario")+ylab("Velocidad (en m/s)")+scale_y_continuous(labels=comma)+ylim(0,15)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 583 rows containing non-finite values (stat_ydensity).

#Variable de dia
BiciMAD$Dia <- weekdays(as.Date(BiciMAD$unplug_hourDate))
BiciMAD$Dia<-as.factor(BiciMAD$Dia)

#Diferencia por dia laborable/fin de semana 
BiciMAD %>% group_by(Dia) %>% summarize(n=n()) %>%  group_by(Dia) %>% summarize(n.m=mean(n)) %>%ggplot(aes(x=Dia,y=n.m)) + geom_bar(stat='identity') + xlab('Dia de la semana') + ylab('Numero de viajes') + ggtitle("Numero de viajes en bicicleta por dia de la semana")+ scale_y_continuous(labels = comma)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)

#Dias respecto a ageRange
ggplot(BiciMAD, aes(x=Dia))+geom_bar(aes(fill=ageRange), position = "fill")+xlab("Dia de la semana")+ylab("Porcentaje(%) de usuarios")

#Dias respecto a user_type
ggplot(BiciMAD, aes(x=Dia))+geom_bar(aes(fill=user_type), position = "fill")+xlab("Dia de la semana")+ylab("Porcentaje(%) de usuarios")

#Dias respecto a travel_time
BiciMAD_paseyocasional<-subset(BiciMAD, user_type==1|user_type==2)

ggplot(BiciMAD_paseyocasional, aes(x=Dia, y=travel_time))+geom_boxplot()+xlab("Dia de la semana")+ylab("Duracion de recorrido (en segundos)")+scale_y_continuous(labels=comma)+ylim(0,2000)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 208077 rows containing non-finite values (stat_boxplot).

rm(BiciMAD_paseyocasional)

#Variable de franja horaria
BiciMAD$unplug_hourTime<-as.numeric(BiciMAD$unplug_hourTime)

BiciMAD$franja_horaria <- ifelse(BiciMAD$unplug_hourTime < 4 | BiciMAD$unplug_hourTime >= 20, "Noche", ifelse(BiciMAD$unplug_hourTime >= 5 & BiciMAD$unplug_hourTime < 12, "Manana", ifelse(BiciMAD$unplug_hourTime >= 12 & BiciMAD$unplug_hourTime < 16, "Tarde", "Tarde-Noche")))

BiciMAD$unplug_hourTime<-as.factor(BiciMAD$unplug_hourTime)
BiciMAD$franja_horaria<-as.factor(BiciMAD$franja_horaria)

#Franja horaria respecto a grupo de edad
ggplot(BiciMAD, aes(x=franja_horaria)) + geom_bar(aes(fill = ageRange), position = "fill")

#Franja horaria respecto a tipo de usuario
ggplot(BiciMAD, aes(x=franja_horaria)) + geom_bar(aes(fill = user_type), position = "fill")

#Variable estaciones
BiciMAD$season <- ifelse( BiciMAD$unplug_hourDate < "2018-03-20" | BiciMAD$unplug_hourDate >="2018-12-21", "Invierno",  ifelse(BiciMAD$unplug_hourDate >="2018-03-20" & BiciMAD$unplug_hourDate < "2018-06-21","Primavera",                  ifelse(BiciMAD$unplug_hourDate >="2018-06-21" & BiciMAD$unplug_hourDate<"2018-09-23","Verano","Otono")))     

BiciMAD$season <-as.factor(BiciMAD$season)

#travel_time de estaciones
ggplot(BiciMAD, aes(season, travel_time)) + geom_boxplot()+xlab("Estacion del ano")+ylab("Duracion del recorrido (en segundos)")+scale_y_continuous(labels = comma)+ylim(0,3000)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 210419 rows containing non-finite values (stat_boxplot).

#20 Marzo a las 16h 14m UTC     -> Primavera  2018-03-20
#21 Junio a las 10h 7m UTC      -> Verano     2018-06-21
#23 Septiembre a las 1h 53m UTC -> Otono      2018-09-23
#21 diciembre a las 22h 22m UTC -> Invierno   2018-12-21

5-Mapas de calor

#Filtraremos unicamente por los usuarios, excluyendo los empleados
BiciMAD_paseyocasional<-subset(BiciMAD, user_type==1|user_type==2)

#Nuevo dataframe, con combinaciones
users_edad_total<-BiciMAD_paseyocasional %>%
  group_by(BiciMAD_paseyocasional$user_type,BiciMAD_paseyocasional$ageRange,BiciMAD_paseyocasional$season,BiciMAD_paseyocasional$unplug_hourDate) %>%
  summarise(Total = n())
## `summarise()` regrouping output by 'BiciMAD_paseyocasional$user_type', 'BiciMAD_paseyocasional$ageRange', 'BiciMAD_paseyocasional$season' (override with `.groups` argument)
#Renombramos
users_edad_total<-setNames(users_edad_total, c("tipo_usuario", "grupo_edad","estacion","Dia", "n"))


#Comparamos grupo de edad con tipo de usuario
ggplot(users_edad_total, aes(x=grupo_edad,y=tipo_usuario))+geom_tile(aes(fill=n))+scale_fill_gradient()

#Comparamos grupo de edad con los dias, para ver quienes son los que mas utilizan las bicicletas a lo largo del ano
ggplot(users_edad_total, aes(x=Dia,y=grupo_edad))+geom_tile(aes(fill=n))+scale_fill_gradient(low="steelBlue",high="white")

rm(users_edad_total)

6-Asimetria por epocas

6.1-Invierno

#Filtramos por la estacion y los usuarios con pase anual y ocasionales, no incluiremos los empleados

BiciMAD_Winter<-subset(BiciMAD, season=="Invierno" & user_type==1 | user_type==2)

head(BiciMAD_Winter[order(BiciMAD_Winter$travel_time, decreasing = T),], n=10)
##               X ageRange idplug_base idplug_station
## 615547   615547        5          20             53
## 217727   217727        3          21             18
## 52137     52137        4          21             79
## 598660   598660        4          12             89
## 286477   286477        5           9             92
## 324594   324594        5          21             29
## 2287342 2287342        0          16            133
## 417552   417552        0          23            117
## 341810   341810        0           1            129
## 3617294 3617294        4           8             46
##                                     address_end  long_end  lat_end
## 615547                    Ronda de Atocha nº 34 -3.699276 40.40609
## 217727          Plaza de Carlos Cambronero nº 2 -3.703831 40.42326
## 52137          Avenida de Menéndez Pelayo nº 90 -3.678484 40.40828
## 598660             Calle Doctor Esquerdo nº 191 -3.672602 40.40325
## 286477                        Calle Ayala nº 44 -3.683257 40.42774
## 324594               Plaza de Celenque nº 1 - B -3.706384 40.41728
## 2287342                Paseo de los Olmos nº 28 -3.710811 40.40341
## 417552       Calle Arcipreste de Hita nº 12 - A -3.717544 40.43373
## 341810  Paseo de Santa María de la Cabeza nº 58 -3.698767 40.40179
## 3617294                      Calle Atocha nº 32 -3.698232 40.41071
##         idunplug_base idunplug_station                            address_start
## 615547             18               91                     Paseo del Prado nº 1
## 217727              7              163 Calle del General Álvarez de Castro nº 2
## 52137              21               91                     Paseo del Prado nº 1
## 598660              2               74         Avenida de Menéndez Pelayo nº 63
## 286477              9               92                        Calle Ayala nº 44
## 324594             21               29               Plaza de Celenque nº 1 - B
## 2287342            10              112                  Calle Serrano nº 34 - B
## 417552             23              117       Calle Arcipreste de Hita nº 12 - A
## 341810             20              129  Paseo de Santa María de la Cabeza nº 58
## 3617294             8               46                       Calle Atocha nº 32
##         long_start lat_start travel_time unplug_hourDate unplug_hourTime
## 615547   -3.692622  40.41861       21595      2018-03-18               6
## 217727   -3.701569  40.43447       21583      2018-01-26              16
## 52137    -3.692622  40.41861       21581      2018-01-09              19
## 598660   -3.677923  40.41531       21572      2018-03-16              13
## 286477   -3.683257  40.42774       21557      2018-02-03               6
## 324594   -3.706384  40.41728       21556      2018-02-09              24
## 2287342  -3.687745  40.42496       21546      2018-09-01               7
## 417552   -3.717544  40.43373       21540      2018-02-19              12
## 341810   -3.698767  40.40179       21538      2018-02-11               7
## 3617294  -3.698232  40.41071       21537      2018-12-22              20
##         user_type zip_code SL_distancia  Velocidad     Dia franja_horaria
## 615547          1    28028     1503.585 0.06962655 domingo         Manana
## 217727          1    28003     1262.336 0.05848752 viernes    Tarde-Noche
## 52137           1    28045     1661.124 0.07697161  martes    Tarde-Noche
## 598660          1    28007     1415.750 0.06562908 viernes          Tarde
## 286477          1    28006        0.000 0.00000000  sábado         Manana
## 324594          1    28109        0.000 0.00000000 viernes          Noche
## 2287342         2     <NA>     3095.106 0.14365105  sábado         Manana
## 417552          1     <NA>        0.000 0.00000000   lunes          Tarde
## 341810          1     <NA>        0.000 0.00000000 domingo         Manana
## 3617294         1     <NA>        0.000 0.00000000  sábado          Noche
##           season
## 615547  Invierno
## 217727  Invierno
## 52137   Invierno
## 598660  Invierno
## 286477  Invierno
## 324594  Invierno
## 2287342   Verano
## 417552  Invierno
## 341810  Invierno
## 3617294 Invierno
NOrigen_winter<-count(BiciMAD_Winter,BiciMAD_Winter$address_start)
NDestino_winter<-count(BiciMAD_Winter,BiciMAD_Winter$address_end)

NOrigen_winter<-setNames(NOrigen_winter, c("Estacion_Origen", "Salidas"))
NDestino_winter<-setNames(NDestino_winter, c("Estacion_Destino", "Entradas"))

Combinado_winter<-left_join(NOrigen_winter, NDestino_winter, by=c("Estacion_Origen"="Estacion_Destino"))

#Creacion de la variable de asimetria
Combinado_winter$Asimetria_total<-Combinado_winter$Entradas - Combinado_winter$Salidas


#Los usuarios anclan en invierno bicis en las siguientes direcciones
head(arrange(Combinado_winter,desc(Asimetria_total)), n = 10)
##                            Estacion_Origen Salidas Entradas Asimetria_total
## 1                Paseo de la Chopera nº 14    7503    10082            2579
## 2                 Paseo de la Florida nº 8    5119     7349            2230
## 3             Calle Cerro de la Plata nº 2    3183     4673            1490
## 4  Paseo de Santa María de la Cabeza nº 58    7815     9054            1239
## 5                 Plaza de la Cebada nº 16    8971     9934             963
## 6                 Paseo de los Olmos nº 28    5172     6071             899
## 7         Calle Jaime el Conquistador nº30    8485     9341             856
## 8           Paseo de las Delicias nº 92-94    4289     4925             636
## 9         Calle Fernando el Católico nº 19    7273     7902             629
## 10        Calle Palos de la Frontera nº 40    5272     5888             616
#Los usuarios desanclan en invierno bicis en las siguientes direcciones
head(arrange(Combinado_winter,Asimetria_total), n = 10)
##                     Estacion_Origen Salidas Entradas Asimetria_total
## 1      Paseo de la Castellana nº 67    3663     2802            -861
## 2       Paseo de la Castellana nº 4    2314     1561            -753
## 3             Calle Carlos III nº 1    4010     3403            -607
## 4       Calle María de Guzmán nº 58    3227     2656            -571
## 5               Calle Valencia nº 1    9439     8905            -534
## 6  Calle Sor Ángela de la Cruz nº 2    2475     1966            -509
## 7              Paseo Yeserías nº 45    2270     1764            -506
## 8   Avenida de Menéndez Pelayo nº 3    3111     2616            -495
## 9           Puerta del Sol nº 1 - B    3386     2893            -493
## 10          Calle Serrano nº 34 - A    1888     1425            -463
#Limpiamos
rm(BiciMAD_Winter)
rm(NOrigen_winter)
rm(NDestino_winter)
rm(Combinado_winter)

6.2-Primavera

#Filtramos por la estacion y los usuarios con pase anual y ocasionales, no incluiremos los empleados

BiciMAD_Spring<-subset(BiciMAD, season=="Primavera" & user_type==1 | user_type==2)

head(BiciMAD_Spring[order(BiciMAD_Spring$travel_time, decreasing = T),], n=10)
##               X ageRange idplug_base idplug_station
## 1484363 1484363        5           3            123
## 1090539 1090539        5          17             17
## 1287925 1287925        3           3            112
## 805356   805356        5           4              8
## 1393845 1393845        5           6             33
## 2287342 2287342        0          16            133
## 1496668 1496668        0           1             17
## 1429500 1429500        5           8            129
## 1168578 1168578        5           2            141
## 718734   718734        4           2             93
##                                     address_end  long_end  lat_end
## 1484363                     Plaza de España - A -3.711603 40.42402
## 1090539                Calle San Bernardo nº 22 -3.707507 40.42307
## 1287925                 Calle Serrano nº 34 - B -3.687745 40.42496
## 805356                    Calle Hortaleza nº 75 -3.697771 40.42519
## 1393845        Calle del Marqués de Cubas nº 25 -3.695736 40.41626
## 2287342                Paseo de los Olmos nº 28 -3.710811 40.40341
## 1496668                Calle San Bernardo nº 22 -3.707507 40.42307
## 1429500 Paseo de Santa María de la Cabeza nº 58 -3.698767 40.40179
## 1168578                      Calle Orense nº 36 -3.694622 40.45485
## 718734                    Calle Velázquez nº 75 -3.683830 40.43136
##         idunplug_base idunplug_station               address_start long_start
## 1484363             3              123         Plaza de España - A  -3.711603
## 1090539            22               55        Ronda de Atocha nº 2  -3.693520
## 1287925            12               10    Plaza de San Miguel nº 9  -3.709508
## 805356              4                8       Calle Hortaleza nº 75  -3.697771
## 1393845             8                1     Puerta del Sol nº 1 - A  -3.702425
## 2287342            10              112     Calle Serrano nº 34 - B  -3.687745
## 1496668             1               17    Calle San Bernardo nº 22  -3.707507
## 1429500            17              109 Paseo de la Castellana nº 4  -3.689534
## 1168578            17              141          Calle Orense nº 36  -3.694622
## 718734              4              160    Calle Cea Bermúdez nº 59  -3.715433
##         lat_start travel_time unplug_hourDate unplug_hourTime user_type
## 1484363  40.42402       21578      2018-06-16               7         1
## 1090539  40.40756       21566      2018-05-08              16         1
## 1287925  40.41561       21558      2018-05-31              24         1
## 805356   40.42519       21555      2018-04-12              22         1
## 1393845  40.41690       21550      2018-06-09              24         1
## 2287342  40.42496       21546      2018-09-01               7         2
## 1496668  40.42307       21538      2018-06-17              12         1
## 1429500  40.42683       21534      2018-06-12               8         1
## 1168578  40.45485       21523      2018-05-15              10         1
## 718734   40.43899       21516      2018-04-02              23         1
##         zip_code SL_distancia  Velocidad     Dia franja_horaria    season
## 1484363    28041       0.0000 0.00000000  sábado         Manana Primavera
## 1090539    28005    2094.4523 0.09711826  martes    Tarde-Noche Primavera
## 1287925    28034    2118.2313 0.09825732  jueves          Noche Primavera
## 805356     28004       0.0000 0.00000000  jueves          Noche Primavera
## 1393845    28031     571.3771 0.02651402  sábado          Noche Primavera
## 2287342     <NA>    3095.1056 0.14365105  sábado         Manana    Verano
## 1496668     <NA>       0.0000 0.00000000 domingo          Tarde Primavera
## 1429500    28005    2895.2560 0.13445045  martes         Manana Primavera
## 1168578    28020       0.0000 0.00000000  martes         Manana Primavera
## 718734     28003    2809.3786 0.13057160   lunes          Noche Primavera
NOrigen_spring<-count(BiciMAD_Spring,BiciMAD_Spring$address_start)
NDestino_spring<-count(BiciMAD_Spring,BiciMAD_Spring$address_end)

NOrigen_spring<-setNames(NOrigen_spring, c("Estacion_Origen", "Salidas"))
NDestino_spring<-setNames(NDestino_spring, c("Estacion_Destino", "Entradas"))

Combinado_spring<-left_join(NOrigen_spring, NDestino_spring, by=c("Estacion_Origen"="Estacion_Destino"))

#Creacion de la variable de asimetria
Combinado_spring$Asimetria_total<-Combinado_spring$Entradas - Combinado_spring$Salidas

#Los usuarios anclan en primavera bicis en las siguientes direcciones
head(arrange(Combinado_spring,desc(Asimetria_total)), n = 10)
##                            Estacion_Origen Salidas Entradas Asimetria_total
## 1                Paseo de la Chopera nº 14   11036    14738            3702
## 2                 Paseo de la Florida nº 8    7206    10074            2868
## 3             Calle Cerro de la Plata nº 2    4679     6528            1849
## 4  Paseo de Santa María de la Cabeza nº 58    9947    11460            1513
## 5                 Paseo de los Olmos nº 28    6746     8074            1328
## 6                 Plaza de la Cebada nº 16   12200    13403            1203
## 7           Paseo de las Delicias nº 92-94    6215     7111             896
## 8         Calle Jaime el Conquistador nº30   11022    11816             794
## 9         Avenida de Menéndez Pelayo nº 38    8501     9201             700
## 10        Calle Fernando el Católico nº 19    9328     9978             650
#Los usuarios desanclan en primavera bicis en las siguientes direcciones
head(arrange(Combinado_spring,Asimetria_total), n = 10)
##                       Estacion_Origen Salidas Entradas Asimetria_total
## 1        Paseo de la Castellana nº 67    4482     3465           -1017
## 2         Paseo de la Castellana nº 4    3209     2265            -944
## 3             Calle Serrano nº 34 - A    2702     1975            -727
## 4     Avenida de Menéndez Pelayo nº 3    4293     3571            -722
## 5               Calle Carlos III nº 1    4824     4133            -691
## 6  Plaza de San Juan de la Cruz nº 11    3380     2699            -681
## 7        Paseo de la Castellana nº 42    4993     4335            -658
## 8                    Calle Goya nº 18    6113     5503            -610
## 9         Calle María de Guzmán nº 58    4222     3627            -595
## 10            Puerta del Sol nº 1 - B    4866     4288            -578
rm(BiciMAD_Spring)
rm(NOrigen_spring)
rm(NDestino_spring)
rm(Combinado_spring)

6.3-Verano

#Filtramos por la estacion y los usuarios con pase anual y ocasionales, no incluiremos los empleados

BiciMAD_Summer<-subset(BiciMAD, season=="Verano" & user_type==1 | user_type==2)

head(BiciMAD_Summer[order(BiciMAD_Summer$travel_time, decreasing = T),], n=10)
##               X ageRange idplug_base idplug_station
## 2243926 2243926        5          19            135
## 2113947 2113947        0          21             99
## 1749101 1749101        0          20            136
## 2094975 2094975        6           1             49
## 2442296 2442296        4          16             87
## 2516073 2516073        0          22             10
## 1935728 1935728        3           4            160
## 1934644 1934644        4          20            163
## 1637717 1637717        5           3             74
## 2349910 2349910        0           1             42
##                                      address_end  long_end  lat_end
## 2243926                Paseo de la Chopera nº 14 -3.697183 40.39194
## 2113947              Paseo de Recoletos nº 32-34 -3.690965 40.42270
## 1749101           Paseo de las Delicias nº 92-94 -3.694503 40.39726
## 2094975     Glorieta de la Puerta de Toledo nº 1 -3.711051 40.40704
## 2442296             Calle Doctor Esquerdo nº 161 -3.669753 40.40846
## 2516073                 Plaza de San Miguel nº 9 -3.709508 40.41561
## 1935728                 Calle Cea Bermúdez nº 59 -3.715433 40.43899
## 1934644 Calle del General Álvarez de Castro nº 2 -3.701569 40.43447
## 1637717         Avenida de Menéndez Pelayo nº 63 -3.677923 40.41531
## 2349910            Carrera de San Francisco nº 1 -3.712073 40.41104
##         idunplug_base idunplug_station                            address_start
## 2243926            19              135                Paseo de la Chopera nº 14
## 2113947            21               99              Paseo de Recoletos nº 32-34
## 1749101            17               81            Plaza del Conde de Casal nº 8
## 2094975             4               49     Glorieta de la Puerta de Toledo nº 1
## 2442296            12               80           Avenida del Mediterráneo nº 19
## 2516073            13              124                      Plaza de España - B
## 1935728            19               18          Plaza de Carlos Cambronero nº 2
## 1934644            20              163 Calle del General Álvarez de Castro nº 2
## 1637717             4               65                    Calle O'Donnell nº 28
## 2349910            21               25                   Calle Jacometrezo nº 3
##         long_start lat_start travel_time unplug_hourDate unplug_hourTime
## 2243926  -3.697183  40.39194       21598      2018-08-25               3
## 2113947  -3.690965  40.42270       21597      2018-08-11              22
## 1749101  -3.670422  40.40635       21595      2018-07-10              12
## 2094975  -3.711051  40.40704       21594      2018-08-09               3
## 2442296  -3.675012  40.40767       21590      2018-09-12              10
## 2516073  -3.711703  40.42412       21576      2018-09-17              10
## 1935728  -3.703831  40.42326       21573      2018-07-25              12
## 1934644  -3.701569  40.43447       21568      2018-07-25              11
## 1637717  -3.675204  40.42140       21559      2018-07-02              24
## 2349910  -3.706538  40.42008       21552      2018-09-06              18
##         user_type zip_code SL_distancia  Velocidad       Dia franja_horaria
## 2243926         1    28011       0.0000 0.00000000    sábado          Noche
## 2113947         1     <NA>       0.0000 0.00000000    sábado          Noche
## 1749101         1     <NA>    2278.3034 0.10550143    martes          Tarde
## 2094975         1    28005       0.0000 0.00000000    jueves          Noche
## 2442296         1    28013     454.2566 0.02104014 miércoles         Manana
## 2516073         1     <NA>     965.8838 0.04476658     lunes         Manana
## 1935728         1    28007    2008.0508 0.09308167 miércoles          Tarde
## 1934644         1    28012       0.0000 0.00000000 miércoles         Manana
## 1637717         1    28032     716.3376 0.03322685     lunes          Noche
## 2349910         1     <NA>    1110.0967 0.05150783    jueves    Tarde-Noche
##         season
## 2243926 Verano
## 2113947 Verano
## 1749101 Verano
## 2094975 Verano
## 2442296 Verano
## 2516073 Verano
## 1935728 Verano
## 1934644 Verano
## 1637717 Verano
## 2349910 Verano
NOrigen_summer<-count(BiciMAD_Summer,BiciMAD_Summer$address_start)
NDestino_summer<-count(BiciMAD_Summer,BiciMAD_Summer$address_end)

NOrigen_summer<-setNames(NOrigen_summer, c("Estacion_Origen", "Salidas"))
NDestino_summer<-setNames(NDestino_summer, c("Estacion_Destino", "Entradas"))

Combinado_summer<-left_join(NOrigen_summer, NDestino_summer, by=c("Estacion_Origen"="Estacion_Destino"))

#Creacion de la variable de asimetria
Combinado_summer$Asimetria_total<-Combinado_summer$Entradas - Combinado_summer$Salidas

#Los usuarios anclan en verano bicis en las siguientes direcciones
head(arrange(Combinado_summer,desc(Asimetria_total)), n = 10)
##                            Estacion_Origen Salidas Entradas Asimetria_total
## 1                Paseo de la Chopera nº 14   12746    17356            4610
## 2                 Paseo de la Florida nº 8    8914    11620            2706
## 3             Calle Cerro de la Plata nº 2    5625     8125            2500
## 4  Paseo de Santa María de la Cabeza nº 58   11057    12910            1853
## 5                 Paseo de los Olmos nº 28    7939     9331            1392
## 6         Calle Jaime el Conquistador nº30   13780    15047            1267
## 7                 Plaza de la Cebada nº 16   12955    14023            1068
## 8           Paseo de las Delicias nº 92-94    7544     8597            1053
## 9         Calle Fernando el Católico nº 19   10366    11060             694
## 10             Plaza de San Francisco nº 5    8111     8801             690
#Los usuarios desanclan en verano bicis en las siguientes direcciones
head(arrange(Combinado_summer,Asimetria_total), n = 10)
##                     Estacion_Origen Salidas Entradas Asimetria_total
## 1      Paseo de la Castellana nº 67    5525     4225           -1300
## 2       Paseo de la Castellana nº 4    3715     2664           -1051
## 3   Avenida de Menéndez Pelayo nº 3    4889     3949            -940
## 4             Calle Carlos III nº 1    5689     4789            -900
## 5           Calle Serrano nº 34 - A    3150     2358            -792
## 6           Puerta del Sol nº 1 - B    5654     4867            -787
## 7           Calle Serrano nº 34 - B    3579     2813            -766
## 8      Paseo de la Castellana nº 42    5614     4890            -724
## 9  Calle Sor Ángela de la Cruz nº 2    4127     3424            -703
## 10              Calle Serrano nº 54    5965     5282            -683
#Limpiamos
rm(BiciMAD_Summer)
rm(NOrigen_summer)
rm(NDestino_summer)
rm(Combinado_summer)

6.4-Otono

#Filtramos por la estacion y los usuarios con pase anual y ocasionales, no incluiremos los empleados

BiciMAD_Fall<-subset(BiciMAD, season=="Otono" & user_type==1 | user_type==2)

head(BiciMAD_Fall[order(BiciMAD_Fall$travel_time, decreasing = T),], n=10)
##               X ageRange idplug_base idplug_station
## 3067369 3067369        0           7            123
## 3293401 3293401        5           5             64
## 2826354 2826354        0           2              5
## 3501301 3501301        4           7             58
## 3214989 3214989        5           1            114
## 3361534 3361534        0          17             31
## 2667824 2667824        5           7             78
## 2812718 2812718        2          20            136
## 2287342 2287342        0          16            133
## 3142926 3142926        4          24             12
##                                 address_end  long_end  lat_end idunplug_base
## 3067369                 Plaza de España - A -3.711603 40.42402             6
## 3293401      Plaza de la Independencia nº 6 -3.688398 40.41975            24
## 2826354         Calle Manuela Malasaña nº 5 -3.702587 40.42855            23
## 3501301                  Calle Barceló nº 7 -3.700423 40.42668            23
## 3214989 Calle Mártires Concepcionistas nº 2 -3.670602 40.42730            18
## 3361534                    Calle Jesús nº 1 -3.695618 40.41328            22
## 2667824    Avenida de Menéndez Pelayo nº 73 -3.676681 40.41176             4
## 2812718      Paseo de las Delicias nº 92-94 -3.694503 40.39726             1
## 2287342            Paseo de los Olmos nº 28 -3.710811 40.40341            10
## 3142926              Calle San Andrés nº 18 -3.703592 40.42695             5
##         idunplug_station                            address_start long_start
## 3067369              149              Calle Santa Engracia nº 168  -3.703667
## 3293401               67                        Calle Ibiza nº 62  -3.670896
## 2826354               67                        Calle Ibiza nº 62  -3.670896
## 3501301              163 Calle del General Álvarez de Castro nº 2  -3.701569
## 3214989              155               Calle María Francisca nº 1  -3.678717
## 3361534              145                       Calle Orense nº 12  -3.695286
## 2667824              155               Calle María Francisca nº 1  -3.678717
## 2812718              114      Calle Mártires Concepcionistas nº 2  -3.670602
## 2287342              112                  Calle Serrano nº 34 - B  -3.687745
## 3142926               40                    Plaza de Ramales nº 1  -3.712404
##         lat_start travel_time unplug_hourDate unplug_hourTime user_type
## 3067369  40.44637       21596      2018-10-27              12         1
## 3293401  40.41792       21574      2018-11-20               8         1
## 2826354  40.41792       21565      2018-10-08              24         1
## 3501301  40.43447       21560      2018-12-11              22         1
## 3214989  40.44423       21556      2018-11-13               4         1
## 3361534  40.44807       21556      2018-11-27              24         1
## 2667824  40.44423       21552      2018-09-27               6         1
## 2812718  40.42730       21552      2018-10-07               1         1
## 2287342  40.42496       21546      2018-09-01               7         2
## 3142926  40.41673       21526      2018-11-05               3         1
##         zip_code SL_distancia  Velocidad     Dia franja_horaria season
## 3067369     <NA>    2576.8885 0.11932249  sábado          Tarde  Otono
## 3293401    28009    1497.2066 0.06939865  martes         Manana  Otono
## 2826354     <NA>    2934.7728 0.13608963   lunes          Noche  Otono
## 3501301    28015     872.6283 0.04047441  martes          Noche  Otono
## 3214989    28002    2005.6407 0.09304327  martes    Tarde-Noche  Otono
## 3361534     <NA>    3872.5064 0.17964865  martes          Noche  Otono
## 2667824    28002    3617.8902 0.16786796  jueves         Manana  Otono
## 2812718    28041    3909.6528 0.18140557 domingo          Noche  Otono
## 2287342     <NA>    3095.1056 0.14365105  sábado         Manana Verano
## 3142926    28014    1360.9072 0.06322156   lunes          Noche  Otono
NOrigen_fall<-count(BiciMAD_Fall,BiciMAD_Fall$address_start)
NDestino_fall<-count(BiciMAD_Fall,BiciMAD_Fall$address_end)

NOrigen_fall<-setNames(NOrigen_fall, c("Estacion_Origen", "Salidas"))
NDestino_fall<-setNames(NDestino_fall, c("Estacion_Destino", "Entradas"))

Combinado_fall<-left_join(NOrigen_fall, NDestino_fall, by=c("Estacion_Origen"="Estacion_Destino"))

#Creacion de la variable de asimetria
Combinado_fall$Asimetria_total<-Combinado_fall$Entradas - Combinado_fall$Salidas


#Los usuarios anclan en otono bicis en las siguientes direcciones
head(arrange(Combinado_fall,desc(Asimetria_total)), n = 15)
##                            Estacion_Origen Salidas Entradas Asimetria_total
## 1                Paseo de la Chopera nº 14   10830    15328            4498
## 2                 Paseo de la Florida nº 8    7377    10100            2723
## 3  Paseo de Santa María de la Cabeza nº 58   10813    13019            2206
## 4             Calle Cerro de la Plata nº 2    5363     7349            1986
## 5         Calle Jaime el Conquistador nº30   12842    14136            1294
## 6                 Plaza de la Cebada nº 16   12287    13489            1202
## 7                 Paseo de los Olmos nº 28    7110     8191            1081
## 8         Calle Fernando el Católico nº 19   10508    11384             876
## 9           Paseo de las Delicias nº 92-94    7100     7880             780
## 10        Avenida de Menéndez Pelayo nº 38    9056     9637             581
## 11   Calle Juan Martín el Empecinado nº 14    4459     5029             570
## 12            Calle Doctor Esquerdo nº 191    3935     4491             556
## 13        Calle Palos de la Frontera nº 40    8569     9104             535
## 14           Carrera de San Francisco nº 1    6110     6605             495
## 15            Avenida de Alfonso XII nº 54    5798     6272             474
#Los usuarios desanclan en otono bicis en las siguientes direcciones
head(arrange(Combinado_fall,Asimetria_total), n = 15)
##                     Estacion_Origen Salidas Entradas Asimetria_total
## 1      Paseo de la Castellana nº 67    4948     3876           -1072
## 2           Calle Serrano nº 34 - A    3082     2080           -1002
## 3   Avenida de Menéndez Pelayo nº 3    4340     3345            -995
## 4       Paseo de la Castellana nº 4    3228     2278            -950
## 5             Calle Carlos III nº 1    5559     4788            -771
## 6           Calle Serrano nº 34 - B    3358     2627            -731
## 7  Calle Sor Ángela de la Cruz nº 2    3958     3234            -724
## 8                  Calle Pavía nº 6    4083     3382            -701
## 9       Calle María de Guzmán nº 58    5154     4463            -691
## 10              Plaza de España - A    4942     4348            -594
## 11      Paseo de Recoletos nº 32-34    5514     4923            -591
## 12             Paseo del Prado nº 1    6425     5848            -577
## 13              Calle Valencia nº 1   13345    12772            -573
## 14                 Calle Goya nº 99    3019     2470            -549
## 15               Calle Orense nº 36    3291     2742            -549
rm(BiciMAD_Fall)
rm(NOrigen_fall)
rm(NDestino_fall)
rm(Combinado_fall)

#Extraemos los resultados y los compararemos.

7-Calculo de asimetria total

#Como podemos observar, los valores mas altos son realizados por usuarios con pase anual, aunque la media de los usuarios ocasionales sea mas alta.

head(BiciMAD[order(BiciMAD$travel_time, decreasing = T),], n=10)
##               X ageRange idplug_base idplug_station
## 2752110 2752110        5           5              4
## 2733274 2733274        0           9            134
## 1653343 1653343        5           8            151
## 122239   122239        1          17            108
## 2827823 2827823        5          13            108
## 3482291 3482291        4          22            154
## 3482290 3482290        5          23            154
## 3482289 3482289        4          24            154
## 2785352 2785352        4          15            111
## 2298733 2298733        0          15            139
##                           address_end  long_end  lat_end idunplug_base
## 2752110 Plaza del Conde Suchil nº 2-4 -3.706917 40.43029            14
## 2733274    Paseo de la Esperanza nº 2 -3.706452 40.40360            17
## 1653343 Paseo de la Castellana nº 106 -3.690861 40.44533             2
## 122239     Calle Ortega y Gasset nº 6 -3.686565 40.43031            15
## 2827823    Calle Ortega y Gasset nº 6 -3.686565 40.43031            24
## 3482291  Paseo de la Castellana nº 67 -3.691793 40.44574            22
## 3482290  Paseo de la Castellana nº 67 -3.691793 40.44574            18
## 3482289  Paseo de la Castellana nº 67 -3.691793 40.44574             6
## 2785352       Calle Serrano nº 34 - A -3.687723 40.42510             7
## 2298733 Calle del General Yagüe nº 57 -3.700967 40.45728            18
##         idunplug_station                     address_start long_start lat_start
## 2752110              131        Calle Guzmán el Bueno nº 2  -3.713341  40.43065
## 2733274               42     Carrera de San Francisco nº 1  -3.712073  40.41104
## 1653343               67                 Calle Ibiza nº 62  -3.670896  40.41792
## 122239               113                Calle Serrano nº 6  -3.688437  40.42152
## 2827823              108        Calle Ortega y Gasset nº 6  -3.686565  40.43031
## 3482291              159 Calle José Gutiérrez Abascal nº 2  -3.690778  40.43968
## 3482290              159 Calle José Gutiérrez Abascal nº 2  -3.690778  40.43968
## 3482289              151     Paseo de la Castellana nº 106  -3.690861  40.44533
## 2785352               91              Paseo del Prado nº 1  -3.692622  40.41861
## 2298733              174              Paseo Yeserías nº 45  -3.708407  40.39741
##         travel_time unplug_hourDate unplug_hourTime user_type zip_code
## 2752110     5730890      2018-10-03              23         3    28003
## 2733274     2833493      2018-10-01              15         3     <NA>
## 1653343     2755282      2018-07-03               5         3     <NA>
## 122239      1651625      2018-01-17               5         3     <NA>
## 2827823     1625256      2018-10-08               5         3     <NA>
## 3482291     1603461      2018-12-09               2         3     <NA>
## 3482290     1603440      2018-12-09               2         3     <NA>
## 3482289     1603113      2018-12-09               2         3     <NA>
## 2785352     1582306      2018-10-05              22         3     <NA>
## 2298733     1346583      2018-09-02              10         3     <NA>
##         SL_distancia    Velocidad       Dia franja_horaria   season
## 2752110    545.75976 9.523124e-05 miércoles          Noche    Otono
## 2733274    955.69720 3.372859e-04     lunes          Tarde    Otono
## 1653343   3488.57197 1.266140e-03    martes         Manana   Verano
## 122239     990.28934 5.995849e-04 miércoles         Manana Invierno
## 2827823      0.00000 0.000000e+00     lunes         Manana    Otono
## 3482291    680.29548 4.242669e-04   domingo          Noche    Otono
## 3482290    680.29548 4.242725e-04   domingo          Noche    Otono
## 3482289     91.25206 5.692179e-05   domingo          Noche    Otono
## 2785352    832.81611 5.263306e-04   viernes          Noche    Otono
## 2298733   6694.50057 4.971473e-03   domingo         Manana   Verano
NOrigen<-count(BiciMAD,BiciMAD$address_start)
NDestino<-count(BiciMAD,BiciMAD$address_end)

NOrigen<-setNames(NOrigen, c("Estacion_Origen", "Salidas"))
NDestino<-setNames(NDestino, c("Estacion_Destino", "Enganchadas"))

Combinado<-left_join(NOrigen, NDestino, by=c("Estacion_Origen"="Estacion_Destino"))

#Creacion de la variable de asimetria
Combinado$Asimetria_total<-Combinado$Enganchadas - Combinado$Salidas

#Crearemos un data frame similar al Combinado con los datos de asimetria, en el cual no se repitan estaciones
datos_estacion<-distinct(BiciMAD, address_start , lat_start , long_start)

#Juntaremos las cifras totales de asimetria con este nuevo dataframe
BiciMAD_Assm<-left_join(Combinado,datos_estacion, by=c("Estacion_Origen"="address_start"))

#Podemos ver cuales son las estaciones con mas desanclajes y anclajes, independientemente de la estacion
head(arrange(BiciMAD_Assm,desc(Asimetria_total)), n = 10)
##                            Estacion_Origen Salidas Enganchadas Asimetria_total
## 1                Paseo de la Chopera nº 14   51926       53274            1348
## 2                 Paseo de la Florida nº 8   34844       35572             728
## 3  Paseo de Santa María de la Cabeza nº 58   48710       49370             660
## 4                 Plaza de la Cebada nº 16   50067       50650             583
## 5             Calle Cerro de la Plata nº 2   26279       26831             552
## 6               Plaza de Jacinto Benavente   30484       30901             417
## 7         Calle Fernando el Católico nº 19   41227       41640             413
## 8         Calle Jaime el Conquistador nº30   50194       50603             409
## 9             Avenida de Alfonso XII nº 54   21408       21753             345
## 10                Paseo de los Olmos nº 28   31976       32302             326
##    long_start lat_start
## 1   -3.697183  40.39194
## 2   -3.721848  40.42208
## 3   -3.698767  40.40179
## 4   -3.708834  40.41127
## 5   -3.674480  40.40091
## 6   -3.703683  40.41468
## 7   -3.708439  40.43385
## 8   -3.698153  40.39620
## 9   -3.688822  40.40981
## 10  -3.710811  40.40341
head(arrange(BiciMAD_Assm,Asimetria_total), n = 10)
##                  Estacion_Origen Salidas Enganchadas Asimetria_total long_start
## 1         Calle Velázquez nº 130   38840       38029            -811  -3.682862
## 2  Calle del General Yagüe nº 57   26722       26044            -678  -3.700967
## 3     Calle Ortega y Gasset nº 6   32224       31623            -601  -3.686565
## 4           Paseo del Prado nº 1   24871       24278            -593  -3.692622
## 5        Calle Serrano nº 34 - A   10961       10382            -579  -3.687723
## 6   Paseo de la Castellana nº 67   18119       17577            -542  -3.691793
## 7    Paseo de la Castellana nº 4   13157       12646            -511  -3.689534
## 8          Calle Velázquez nº 75   18376       17876            -500  -3.683830
## 9        Calle Serrano nº 34 - B   12501       12087            -414  -3.687745
## 10              Calle Goya nº 18   25666       25264            -402  -3.683788
##    lat_start
## 1   40.43794
## 2   40.45728
## 3   40.43031
## 4   40.41861
## 5   40.42510
## 6   40.44574
## 7   40.42683
## 8   40.43136
## 9   40.42496
## 10  40.42504
#Limpiamos
rm(NDestino)
rm(NOrigen)
rm(datos_estacion)
rm(Combinado)

BiciMAD$X<-NULL

8-Rutas mas populares

8.1-Rutas mas populares por usuarios suscriptores y ocasionales

#Para determinar el comportamiento de los usuarios realizaremos un estudio conjunto de los usuarios user_type = 1 y user_type = 2, aunque realizaremos distinciones entre ambos grupos si es necesario


BiciMAD_paseyocasional$ruta <- paste(BiciMAD_paseyocasional$address_start, BiciMAD_paseyocasional$address_end, sep= " -> ")

#Rutas mas populares por los usuarios
BiciMAD_paseyocasional%>% group_by(ruta) %>% summarize(count=n()) %>% arrange(desc(count)) %>% top_n(n=20) %>% knitr::kable(format="markdown")
## `summarise()` ungrouping output (override with `.groups` argument)
## Selecting by count
ruta count
Plaza de la Independencia nº 6 -> Plaza de la Independencia nº 6 3001
Plaza de Alonso Martínez nº 5 -> Calle Santa Engracia nº 168 2396
Paseo de la Chopera nº 14 -> Paseo de la Chopera nº 14 2393
Paseo de la Florida nº 8 -> Paseo de la Florida nº 8 2127
Calle Valencia nº 1 -> Paseo de la Chopera nº 14 1950
Paseo de la Chopera nº 14 -> Paseo de la Florida nº 8 1908
Paseo de la Florida nº 8 -> Paseo de la Chopera nº 14 1766
Avenida de Menéndez Pelayo nº 11 -> Avenida de Menéndez Pelayo nº 11 1740
Glorieta de Embajadores nº 2 -> Paseo de la Chopera nº 14 1707
Calle Santa Engracia nº 14 -> Calle Santa Engracia nº 168 1667
Calle Jaime el Conquistador nº30 -> Plaza de la Cebada nº 16 1652
Calle Jaime el Conquistador nº30 -> Paseo de la Chopera nº 14 1646
Calle Barceló nº 7 -> Calle Barceló nº 7 1642
Calle del General Álvarez de Castro nº 2 -> Calle Santa Engracia nº 168 1616
Paseo de la Chopera nº 14 -> Calle Jaime el Conquistador nº30 1578
Calle del General Álvarez de Castro nº 2 -> Calle Cea Bermúdez nº 59 1565
Calle Valencia nº 1 -> Calle Valencia nº 1 1522
Calle Cea Bermúdez nº 59 -> Calle del General Álvarez de Castro nº 2 1503
Glorieta de la Puerta de Toledo nº 1 -> Paseo de la Chopera nº 14 1503
Plaza de la Cebada nº 16 -> Paseo de la Chopera nº 14 1493 
#Estaciones de inicio mas populares por parte de los usuarios 
BiciMAD_paseyocasional%>% group_by(address_start) %>% summarize(count=n()) %>% arrange(desc(count)) %>% top_n(n=20) %>% knitr::kable(format="markdown")
## `summarise()` ungrouping output (override with `.groups` argument)
## Selecting by count
address_start count
Calle Valencia nº 1 47055
Calle del General Álvarez de Castro nº 2 46480
Calle Jaime el Conquistador nº30 44320
Plaza de la Cebada nº 16 44091
Calle Santa Engracia nº 168 41723
Avenida de Menéndez Pelayo nº 11 40612
Paseo de Santa María de la Cabeza nº 58 39008
Paseo de la Chopera nº 14 38290
Calle Cea Bermúdez nº 59 36761
Calle Fernando el Católico nº 19 36611
Calle Barceló nº 7 36481
Plaza de la Independencia nº 6 34486
Plaza de Alonso Martínez nº 5 33528
Plaza de Pedro Zerolo 33181
Avenida de Menéndez Pelayo nº 38 33027
Puerta del Sol nº 1 - A 32317
Glorieta de la Puerta de Toledo nº 1 32052
Calle Velázquez nº 130 31600
Calle Juan Bravo nº 50 31018
Plaza de San Ildefonso nº 3 30995 
#Estaciones de destino mas populares por parte de los usuarios
BiciMAD_paseyocasional%>% group_by(address_end) %>% summarize(count=n()) %>% arrange(desc(count)) %>% top_n(n=20) %>% knitr::kable(format="markdown")
## `summarise()` ungrouping output (override with `.groups` argument)
## Selecting by count
address_end count
Paseo de la Chopera nº 14 52677
Calle Jaime el Conquistador nº30 48471
Plaza de la Cebada nº 16 48101
Calle del General Álvarez de Castro nº 2 47203
Paseo de Santa María de la Cabeza nº 58 45738
Calle Valencia nº 1 45051
Calle Santa Engracia nº 168 41780
Avenida de Menéndez Pelayo nº 11 39520
Calle Fernando el Católico nº 19 39430
Calle Cea Bermúdez nº 59 37340
Plaza de la Independencia nº 6 35818
Calle Barceló nº 7 35502
Avenida de Menéndez Pelayo nº 38 35486
Paseo de la Florida nº 8 35189
Glorieta de la Puerta de Toledo nº 1 33577
Plaza de Alonso Martínez nº 5 33278
Plaza de Pedro Zerolo 33119
Calle Velázquez nº 130 32303
Calle Palos de la Frontera nº 40 31548
Plaza de San Ildefonso nº 3 31162

8.2-Rutas mas populares por los empleados

#Para determinar el comportamiento de los empleados realizaremos un estudio con los datos user_type = 3

BiciMAD_empleados<-subset(BiciMAD, user_type==3 & address_end!=address_start)

BiciMAD_empleados$ruta <- paste(BiciMAD_empleados$address_start, BiciMAD_empleados$address_end, sep= " -> ")


#Rutas mas populares por los empleados
BiciMAD_empleados%>% group_by(ruta) %>% summarize(count=n()) %>% arrange(desc(count)) %>% top_n(n=20) %>% knitr::kable(format="markdown")
## `summarise()` ungrouping output (override with `.groups` argument)
## Selecting by count
ruta count
Paseo de la Chopera nº 14 -> Calle Velázquez nº 130 555
Paseo de Santa María de la Cabeza nº 58 -> Calle Velázquez nº 130 484
Paseo de la Chopera nº 14 -> Paseo del Prado nº 1 412
Paseo de la Chopera nº 14 -> Paseo de la Castellana nº 42 404
Paseo de la Chopera nº 14 -> Paseo de la Castellana nº 4 398
Paseo de Santa María de la Cabeza nº 58 -> Calle Goya nº 18 394
Paseo de la Chopera nº 14 -> Paseo de Recoletos nº 32-34 365
Calle Cerro de la Plata nº 2 -> Calle Velázquez nº 130 349
Paseo de la Chopera nº 14 -> Calle Goya nº 18 344
Paseo de la Chopera nº 14 -> Calle Velázquez nº 75 341
Paseo de la Chopera nº 14 -> Plaza de la Independencia nº 6 341
Paseo de Santa María de la Cabeza nº 58 -> Calle Serrano nº 34 - A 340
Paseo de la Florida nº 8 -> Calle Goya nº 18 339
Paseo de Santa María de la Cabeza nº 58 -> Paseo de Recoletos nº 32-34 338
Calle del General Yagüe nº 57 -> Calle Orense nº 36 334
Plaza de la Independencia nº 6 -> Avenida de Menéndez Pelayo nº 90 331
Calle Santa Engracia nº 168 -> Paseo de la Castellana nº 67 325
Paseo de la Chopera nº 14 -> Calle Ortega y Gasset nº 6 323
Calle Cerro de la Plata nº 2 -> Calle Goya nº 18 321
Paseo de Santa María de la Cabeza nº 58 -> Paseo de la Castellana nº 4 320 
#Estaciones de inicio mas populares por parte de los empleados 
BiciMAD_empleados%>% group_by(address_start) %>% summarize(count=n()) %>% arrange(desc(count)) %>% top_n(n=20) %>% knitr::kable(format="markdown")
## `summarise()` ungrouping output (override with `.groups` argument)
## Selecting by count
address_start count
Paseo de la Chopera nº 14 13452
Paseo de Santa María de la Cabeza nº 58 9636
Paseo de la Florida nº 8 9364
Calle Cerro de la Plata nº 2 7643
Calle Velázquez nº 130 7116
Paseo de los Olmos nº 28 6101
Plaza de la Cebada nº 16 5910
Calle Jaime el Conquistador nº30 5800
Avenida de Menéndez Pelayo nº 38 5254
Plaza de la Independencia nº 6 5226
Avenida de Menéndez Pelayo nº 90 4695
Calle Fernando el Católico nº 19 4532
Calle Ortega y Gasset nº 6 4409
Paseo de las Delicias nº 92-94 4267
Calle Palos de la Frontera nº 40 3950
Carrera de San Francisco nº 1 3766
Plaza de San Francisco nº 5 3528
Paseo del Prado nº 1 3511
Calle Antonio Maura nº 15 3368
Paseo de la Esperanza nº 2 3182 
#Estaciones de destino mas populares por parte de los empleados 
BiciMAD_empleados%>% group_by(address_end) %>% summarize(count=n()) %>% arrange(desc(count)) %>% top_n(n=20) %>% knitr::kable(format="markdown")
## `summarise()` ungrouping output (override with `.groups` argument)
## Selecting by count
address_end count
Calle Velázquez nº 130 5630
Paseo de la Castellana nº 4 4346
Calle Ortega y Gasset nº 6 4328
Calle Goya nº 18 3995
Paseo de la Castellana nº 67 3842
Paseo de la Castellana nº 42 3669
Paseo del Prado nº 1 3620
Plaza de la Independencia nº 6 3598
Paseo de Santa María de la Cabeza nº 58 3545
Calle Velázquez nº 75 3286
Avenida de Menéndez Pelayo nº 3 3171
Avenida de Menéndez Pelayo nº 90 3062
Calle Serrano nº 6 3030
Avenida de Menéndez Pelayo nº 38 2923
Calle Serrano nº 34 - A 2915
Calle Valencia nº 1 2900
Calle Serrano nº 54 2863
Paseo de Recoletos nº 32-34 2849
Calle Antonio Maura nº 15 2776
Puerta del Sol nº 1 - A 2768

9-Mapas de Densidad y Asimetria

9.1-Mapa: Densidad de bicicletas Salidas

#register_google(key = "Tu API de Google Cloud", write=TRUE)

#25% como muestra sin utilizar el paquete caret
tamano_muestra <- floor(0.25 * nrow(BiciMAD))

#Metemos la semilla 1 para reproducir la misma reasignacion en otros ordenadores
set.seed(1)
index_muestra<- sample(seq_len(nrow(BiciMAD)), size = tamano_muestra)

muestra <- BiciMAD[index_muestra, ]

rm(train_ind)
## Warning in rm(train_ind): objeto 'train_ind' no encontrado
rm(index_muestra)
rm(tamano_muestra)



BiciMAD_fill<-BiciMAD_paseyocasional %>% group_by(long_start, lat_start) %>% summarise(Total=n())
## `summarise()` regrouping output by 'long_start' (override with `.groups` argument)
#Sacamos el mapa de Madrid de Google Maps mediante nuestra API y ggmap
mapa_densidad<-ggmap(get_googlemap(center = c(lon = -3.6859308, lat = 40.4428951), zoom = 13, scale = 2, maptype ="terrain", color = "color")) 
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.442895,-3.685931&zoom=13&size=640x640&scale=2&maptype=terrain&key=xxx
#Areas de madrid con mas desenganches sobre una muestra del 25% del dataset
mapa_densidad + stat_density2d(aes(x = long_start , y = lat_start, fill = ..level.., alpha =0.2), data = muestra, geom = "polygon") 
## Warning: Removed 46317 rows containing non-finite values (stat_density2d).

mapa_densidad + stat_density2d(aes(x = long_start , y = lat_start, fill = ..level.., alpha =0.1), data = BiciMAD_fill, geom = "polygon") 
## Warning: Removed 6 rows containing non-finite values (stat_density2d).

rm(muestra)
rm(BiciMAD_fill)

9.2-Mapa: Asimetria total

BiciMAD_Assm$long_start<-as.numeric(as.character(BiciMAD_Assm$long_start))
BiciMAD_Assm$lat_start<-as.numeric(as.character(BiciMAD_Assm$lat_start))

#Creacion de la variable de asimetria
BiciMAD_Assm$Categoria_Asimetria<-ifelse(BiciMAD_Assm$Asimetria_total<=-800, "Muy Critico (< -800)",ifelse(BiciMAD_Assm$Asimetria_total >-800 & BiciMAD_Assm$Asimetria_total<=-500, "Critico (-800 to -500)",ifelse(BiciMAD_Assm$Asimetria_total>-500 & BiciMAD_Assm$Asimetria_total<=0, "Algo Critico(-500 to 0)",ifelse(BiciMAD_Assm$Asimetria_total>=0 & BiciMAD_Assm$Asimetria_total<10, "Normal(10-0)","Buena(>10)"))))

#funcion de paleta de colores
ColorPalette<- colorFactor(palette=c("Red","Blue","Orange", "Black", "Green"),domain= BiciMAD_Assm$Categoria_Asimetria)

#Mapa de asimetria
BiciMAD_mapa <-leaflet(data=BiciMAD_Assm) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  setView(lng=-3.6859308,lat=40.4428951, zoom = 12) %>% 
  addCircleMarkers(lng=BiciMAD_Assm$long_start,lat=BiciMAD_Assm$lat_start,
                   clusterOptions = markerClusterOptions(),         
                   color=ColorPalette(BiciMAD_Assm$Categoria_Asimetria),
                   label=paste("Nombre:",BiciMAD_Assm$Estacion_Origen,". Asim:", BiciMAD_Assm$Asimetria_total))%>%
  addLegend(position="bottomright", pal=ColorPalette, values=BiciMAD_Assm$Categoria_Asimetria, title="Nivel de asimetria")

#Visualizacion del mapa
BiciMAD_mapa

9.3-Mapa: Top 10 estaciones con mas y menos asimetria total

#Top 10 estaciones con asimetria negativa
top10negativo<-head(arrange(BiciMAD_Assm,Asimetria_total), n = 10)
top10negativo
##                  Estacion_Origen Salidas Enganchadas Asimetria_total long_start
## 1         Calle Velázquez nº 130   38840       38029            -811  -3.682862
## 2  Calle del General Yagüe nº 57   26722       26044            -678  -3.700967
## 3     Calle Ortega y Gasset nº 6   32224       31623            -601  -3.686565
## 4           Paseo del Prado nº 1   24871       24278            -593  -3.692622
## 5        Calle Serrano nº 34 - A   10961       10382            -579  -3.687723
## 6   Paseo de la Castellana nº 67   18119       17577            -542  -3.691793
## 7    Paseo de la Castellana nº 4   13157       12646            -511  -3.689534
## 8          Calle Velázquez nº 75   18376       17876            -500  -3.683830
## 9        Calle Serrano nº 34 - B   12501       12087            -414  -3.687745
## 10              Calle Goya nº 18   25666       25264            -402  -3.683788
##    lat_start     Categoria_Asimetria
## 1   40.43794    Muy Critico (< -800)
## 2   40.45728  Critico (-800 to -500)
## 3   40.43031  Critico (-800 to -500)
## 4   40.41861  Critico (-800 to -500)
## 5   40.42510  Critico (-800 to -500)
## 6   40.44574  Critico (-800 to -500)
## 7   40.42683  Critico (-800 to -500)
## 8   40.43136  Critico (-800 to -500)
## 9   40.42496 Algo Critico(-500 to 0)
## 10  40.42504 Algo Critico(-500 to 0)
#Top 10 estaciones con asimetria positiva
top10positivo<-head(arrange(BiciMAD_Assm,desc(Asimetria_total)), n = 10)
top10positivo
##                            Estacion_Origen Salidas Enganchadas Asimetria_total
## 1                Paseo de la Chopera nº 14   51926       53274            1348
## 2                 Paseo de la Florida nº 8   34844       35572             728
## 3  Paseo de Santa María de la Cabeza nº 58   48710       49370             660
## 4                 Plaza de la Cebada nº 16   50067       50650             583
## 5             Calle Cerro de la Plata nº 2   26279       26831             552
## 6               Plaza de Jacinto Benavente   30484       30901             417
## 7         Calle Fernando el Católico nº 19   41227       41640             413
## 8         Calle Jaime el Conquistador nº30   50194       50603             409
## 9             Avenida de Alfonso XII nº 54   21408       21753             345
## 10                Paseo de los Olmos nº 28   31976       32302             326
##    long_start lat_start Categoria_Asimetria
## 1   -3.697183  40.39194          Buena(>10)
## 2   -3.721848  40.42208          Buena(>10)
## 3   -3.698767  40.40179          Buena(>10)
## 4   -3.708834  40.41127          Buena(>10)
## 5   -3.674480  40.40091          Buena(>10)
## 6   -3.703683  40.41468          Buena(>10)
## 7   -3.708439  40.43385          Buena(>10)
## 8   -3.698153  40.39620          Buena(>10)
## 9   -3.688822  40.40981          Buena(>10)
## 10  -3.710811  40.40341          Buena(>10)
#Las juntamos en un dataframe
top10combinado<-full_join(top10negativo,top10positivo)
## Joining, by = c("Estacion_Origen", "Salidas", "Enganchadas", "Asimetria_total", "long_start", "lat_start", "Categoria_Asimetria")
top10combinado
##                            Estacion_Origen Salidas Enganchadas Asimetria_total
## 1                   Calle Velázquez nº 130   38840       38029            -811
## 2            Calle del General Yagüe nº 57   26722       26044            -678
## 3               Calle Ortega y Gasset nº 6   32224       31623            -601
## 4                     Paseo del Prado nº 1   24871       24278            -593
## 5                  Calle Serrano nº 34 - A   10961       10382            -579
## 6             Paseo de la Castellana nº 67   18119       17577            -542
## 7              Paseo de la Castellana nº 4   13157       12646            -511
## 8                    Calle Velázquez nº 75   18376       17876            -500
## 9                  Calle Serrano nº 34 - B   12501       12087            -414
## 10                        Calle Goya nº 18   25666       25264            -402
## 11               Paseo de la Chopera nº 14   51926       53274            1348
## 12                Paseo de la Florida nº 8   34844       35572             728
## 13 Paseo de Santa María de la Cabeza nº 58   48710       49370             660
## 14                Plaza de la Cebada nº 16   50067       50650             583
## 15            Calle Cerro de la Plata nº 2   26279       26831             552
## 16              Plaza de Jacinto Benavente   30484       30901             417
## 17        Calle Fernando el Católico nº 19   41227       41640             413
## 18        Calle Jaime el Conquistador nº30   50194       50603             409
## 19            Avenida de Alfonso XII nº 54   21408       21753             345
## 20                Paseo de los Olmos nº 28   31976       32302             326
##    long_start lat_start     Categoria_Asimetria
## 1   -3.682862  40.43794    Muy Critico (< -800)
## 2   -3.700967  40.45728  Critico (-800 to -500)
## 3   -3.686565  40.43031  Critico (-800 to -500)
## 4   -3.692622  40.41861  Critico (-800 to -500)
## 5   -3.687723  40.42510  Critico (-800 to -500)
## 6   -3.691793  40.44574  Critico (-800 to -500)
## 7   -3.689534  40.42683  Critico (-800 to -500)
## 8   -3.683830  40.43136  Critico (-800 to -500)
## 9   -3.687745  40.42496 Algo Critico(-500 to 0)
## 10  -3.683788  40.42504 Algo Critico(-500 to 0)
## 11  -3.697183  40.39194              Buena(>10)
## 12  -3.721848  40.42208              Buena(>10)
## 13  -3.698767  40.40179              Buena(>10)
## 14  -3.708834  40.41127              Buena(>10)
## 15  -3.674480  40.40091              Buena(>10)
## 16  -3.703683  40.41468              Buena(>10)
## 17  -3.708439  40.43385              Buena(>10)
## 18  -3.698153  40.39620              Buena(>10)
## 19  -3.688822  40.40981              Buena(>10)
## 20  -3.710811  40.40341              Buena(>10)
#Creamos una nueva variable de categoria de asimetria
top10combinado$Categoria_Asimetria<-ifelse(top10combinado$Asimetria_total<0,"Asimetria Negativa","Asimetria Positiva")

#Creamos una nueva paleta de colores
ColorPaleta_top10<- colorFactor(palette=c("Red","Green"),domain= top10combinado$Categoria_Asimetria)

#Construiremos el mapa
Top10_mapa <- leaflet(data=top10combinado) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%     
  setView(lng=-3.6859308,lat=40.4428951, zoom = 12) %>%                           
  addCircleMarkers(lng=top10combinado$long_start,lat=top10combinado$lat_start,
                   clusterOptions = markerClusterOptions(),
                   color=ColorPaleta_top10(top10combinado$Categoria_Asimetria),
                   label=paste("Nombre:",top10combinado$Estacion_Origen,". Asim:", top10combinado$Asimetria_total))%>%
  addLegend(position="bottomright", pal=ColorPaleta_top10, values=top10combinado$Categoria_Asimetria, title="Nivel de asimetria")

#Visualizamos
Top10_mapa
#utils::write.csv(BiciMAD_Assm,"BiciMAD_Assm.csv")
#utils::write.csv(BiciMAD,"BiciMAD_clean.csv")
#utils::write.csv(top10negativo,"top10negativo.csv")
#utils::write.csv(top10positivo,"top10positivo.csv")

rm(BiciMAD_Assm)
rm(BiciMAD)
rm(mapa_densidad)
rm(ColorPalette)
rm(BiciMAD_mapa)
rm(ColorPaleta_top10)
rm(Top10_mapa)
rm(top10negativo)
rm(top10positivo)

10-Regresion lineal: Demanda vs.clima

10.1-Agrupacion del numero de usuarios respecto a las variables elegidas

#Contamos el numero de bicicletas por dia/demanda diaria

Bicis_dia<-count(BiciMAD_paseyocasional,unplug_hourDate)

Bicis_dia<-setNames(Bicis_dia,c("Date","Demanda"))
Bicis_dia$Date<-as.Date(Bicis_dia$Date,format="%Y/%m/%d")

#Leemos los datos del clima para 2018. Utilizamos la estacion del Parque de El Retiro ya que esta proxima a muchas estaciones de BiciMAD, es una de las zonas con mas estaciones proximas y la estacion meteorologica que mejor mide, ya que esta poco expuesta a alteraciones por radiacion, debido a encontrarse en una zona verde.
datos_clima<-read.csv("2018 Weatherdata_RETIRO.csv", sep=";")
datos_clima$Date<-as.Date(datos_clima$Date,format="%d/%m/%Y")


#Juntaremos los datos
datos_modelo<-merge(Bicis_dia, datos_clima, by="Date", all.x = TRUE, all.y=TRUE)

#Quitamos los NAs en la variable demanda, debido a que BiciMAD no esta disponible todos los dias del ano
datos_modelo<-datos_modelo[!is.na(datos_modelo$Demanda), ]

#Quitaremos las variables irrelevantes
datos_modelo$Station_id<-NULL
datos_modelo$Station_name<-NULL
datos_modelo$Province<-NULL
datos_modelo$Height_in_meters<-NULL


#Formato variables a numericas
datos_modelo$Rain_amount<-as.numeric(as.character(datos_modelo$Rain_amount))
## Warning: NAs introducidos por coerción
datos_modelo$Tmean<-as.numeric(as.character(datos_modelo$Tmean))
datos_modelo$Tmin<-as.numeric(as.character(datos_modelo$Tmin))
datos_modelo$Wind_Direction<-as.numeric(as.character(datos_modelo$Wind_Direction))
datos_modelo$Wind_speed<-as.numeric(as.character(datos_modelo$Wind_speed))
datos_modelo$Wind_Gust<-as.numeric(as.character(datos_modelo$Wind_Gust))
datos_modelo$Pressure_max<-as.numeric(as.character(datos_modelo$Pressure_max))
datos_modelo$Pressure_min<-as.numeric(as.character(datos_modelo$Pressure_min))
datos_modelo$Month<-as.numeric(as.character(datos_modelo$Month))
datos_modelo$Tmax<-as.numeric(as.character(datos_modelo$Tmax))

10.2-Regresiones lineales simples

A simple vista, no parece que una regresion lineal con relaciones no lineales sea lo ideal. Decidimos dibujar el siguiente modelo lineal simple y construimos uno mejor en base a esto

library(stats)
summary(datos_modelo)
##       Date               Demanda           Day            Month       
##  Min.   :2018-01-01   Min.   :    8   Min.   : 1.00   Min.   : 1.000  
##  1st Qu.:2018-03-31   1st Qu.: 7376   1st Qu.: 8.00   1st Qu.: 3.000  
##  Median :2018-07-05   Median : 9351   Median :15.00   Median : 7.000  
##  Mean   :2018-07-02   Mean   : 9490   Mean   :15.35   Mean   : 6.524  
##  3rd Qu.:2018-10-04   3rd Qu.:11913   3rd Qu.:23.00   3rd Qu.:10.000  
##  Max.   :2019-01-01   Max.   :16562   Max.   :31.00   Max.   :12.000  
##                                                                       
##       Year          Tmean        Rain_amount          Tmin         Time_tmin  
##  Min.   :2018   Min.   : 1.60   Min.   : 0.000   Min.   :-1.90   Varias : 30  
##  1st Qu.:2018   1st Qu.: 8.20   1st Qu.: 0.000   1st Qu.: 4.50   23:59  : 27  
##  Median :2018   Median :14.20   Median : 0.000   Median :10.00   6:30   : 24  
##  Mean   :2018   Mean   :15.38   Mean   : 1.619   Mean   :10.63   5:30   : 22  
##  3rd Qu.:2018   3rd Qu.:23.00   3rd Qu.: 0.700   3rd Qu.:17.20   6:10   : 16  
##  Max.   :2019   Max.   :32.90   Max.   :28.900   Max.   :25.90   5:10   : 14  
##                                 NA's   :5                        (Other):224  
##       Tmax         Time_tmax   Wind_Direction    Wind_speed      Wind_Gust     
##  Min.   : 2.80   13:50  : 29   Min.   : 1.00   Min.   :0.000   Min.   : 2.800  
##  1st Qu.:12.00   13:30  : 28   1st Qu.: 5.00   1st Qu.:1.100   1st Qu.: 6.900  
##  Median :18.50   14:30  : 23   Median :19.00   Median :1.700   Median : 8.900  
##  Mean   :20.13   14:40  : 21   Mean   :20.23   Mean   :1.933   Mean   : 9.119  
##  3rd Qu.:28.50   14:10  : 20   3rd Qu.:26.00   3rd Qu.:2.500   3rd Qu.:11.100  
##  Max.   :39.90   14:50  : 16   Max.   :99.00   Max.   :6.100   Max.   :23.300  
##                  (Other):220   NA's   :18      NA's   :17      NA's   :18      
##    Time_Gust    Pressure_max   Time_pressure_max  Pressure_min  
##         : 18   Min.   :924.4   24     :69        Min.   :915.0  
##  Varias : 16   1st Qu.:939.6   0      :63        1st Qu.:935.7  
##  0:10   :  9   Median :942.1   Varias :43        Median :938.2  
##  13:30  :  8   Mean   :942.6   10     :42        Mean   :938.3  
##  16:50  :  8   3rd Qu.:946.4          :31        3rd Qu.:941.7  
##  0:00   :  7   Max.   :961.9   9      :22        Max.   :953.7  
##  (Other):291   NA's   :31      (Other):87        NA's   :31     
##  Time_pressure_min
##  Varias : 51      
##  17     : 41      
##  18     : 41      
##  0      : 35      
##  24     : 34      
##         : 31      
##  (Other):124
#Primer modelo linear
modelo<-lm(Demanda~Month+Tmean+Tmin+Tmax+Rain_amount+Wind_Direction+Wind_Gust+Wind_speed+Pressure_min+Pressure_max,data=datos_modelo)
summary(modelo)
## 
## Call:
## lm(formula = Demanda ~ Month + Tmean + Tmin + Tmax + Rain_amount + 
##     Wind_Direction + Wind_Gust + Wind_speed + Pressure_min + 
##     Pressure_max, data = datos_modelo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10800.9  -1802.3    393.7   1775.9   4704.4 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    38726.289  30553.189   1.268 0.205982    
## Month             90.025     49.763   1.809 0.071466 .  
## Tmean          -2210.548   4513.489  -0.490 0.624667    
## Tmin            1124.228   2262.465   0.497 0.619629    
## Tmax            1238.300   2255.783   0.549 0.583462    
## Rain_amount     -197.977     59.140  -3.348 0.000922 ***
## Wind_Direction     1.208      7.738   0.156 0.876021    
## Wind_Gust       -156.045     84.628  -1.844 0.066209 .  
## Wind_speed        75.515    249.823   0.302 0.762658    
## Pressure_min     -25.293     74.802  -0.338 0.735502    
## Pressure_max      -7.734     77.260  -0.100 0.920336    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2628 on 293 degrees of freedom
##   (53 observations deleted due to missingness)
## Multiple R-squared:  0.3162, Adjusted R-squared:  0.2929 
## F-statistic: 13.55 on 10 and 293 DF,  p-value: < 2.2e-16
#Segundo modelo linear
modelo2<-lm(Demanda~Month+Tmean+Rain_amount+Wind_Gust+Wind_speed+Pressure_min+Pressure_max,data=datos_modelo)
summary(modelo2)
## 
## Call:
## lm(formula = Demanda ~ Month + Tmean + Rain_amount + Wind_Gust + 
##     Wind_speed + Pressure_min + Pressure_max, data = datos_modelo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10993.7  -1760.7    375.7   1776.3   4794.6 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39165.632  30116.201   1.300 0.194447    
## Month           77.237     46.106   1.675 0.094948 .  
## Tmean          165.597     21.419   7.731 1.67e-13 ***
## Rain_amount   -208.253     56.568  -3.681 0.000276 ***
## Wind_Gust     -153.952     83.650  -1.840 0.066706 .  
## Wind_speed      47.135    241.647   0.195 0.845483    
## Pressure_min   -25.802     74.127  -0.348 0.728033    
## Pressure_max    -7.179     76.836  -0.093 0.925622    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2619 on 296 degrees of freedom
##   (53 observations deleted due to missingness)
## Multiple R-squared:  0.3142, Adjusted R-squared:  0.2979 
## F-statistic: 19.37 on 7 and 296 DF,  p-value: < 2.2e-16
#Tercer modelo linear
modelo3<-lm(Demanda~Month+Tmean+Rain_amount+Wind_Gust,data=datos_modelo)
summary(modelo3)
## 
## Call:
## lm(formula = Demanda ~ Month + Tmean + Rain_amount + Wind_Gust, 
##     data = datos_modelo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11019.7  -1685.4    339.8   1865.2   4481.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8030.63     632.16  12.703  < 2e-16 ***
## Month          85.97      43.43   1.980  0.04859 *  
## Tmean         174.37      18.61   9.370  < 2e-16 ***
## Rain_amount  -181.66      39.59  -4.588 6.37e-06 ***
## Wind_Gust    -162.20      50.95  -3.183  0.00159 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2548 on 329 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.3825, Adjusted R-squared:  0.375 
## F-statistic: 50.95 on 4 and 329 DF,  p-value: < 2.2e-16
#Introducimos ruido a los datos con NA
datos_modelo$Wind_Gust[is.na(datos_modelo$Wind_Gust)] <- 0 
datos_modelo$Rain_amount[is.na(datos_modelo$Rain_amount)] <- 0 

#Este modelo sera peor que el 5 debido al ruido, pero podremos extraer los resultados al tener los datos del modelo y el modelo6, el mismo numero de registros/observaciones en total
modelo6<-lm(Demanda~Month+Tmean+Rain_amount+Wind_Gust,data=datos_modelo)
summary(modelo6)
## 
## Call:
## lm(formula = Demanda ~ Month + Tmean + Rain_amount + Wind_Gust, 
##     data = datos_modelo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11251.5  -1609.0    444.1   1831.4   6906.4 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7104.29     541.69  13.115  < 2e-16 ***
## Month         114.49      42.97   2.664  0.00807 ** 
## Tmean         173.18      18.74   9.240  < 2e-16 ***
## Rain_amount  -179.24      34.92  -5.132 4.74e-07 ***
## Wind_Gust     -85.33      40.22  -2.122  0.03455 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2652 on 352 degrees of freedom
## Multiple R-squared:  0.3486, Adjusted R-squared:  0.3412 
## F-statistic:  47.1 on 4 and 352 DF,  p-value: < 2.2e-16
#Extraeremos el modelo final 
Fitted<-predict(modelo6)
Fitted<-as.data.frame(Fitted)
datos_modelo$Fitted<- Fitted$Fitted


#Con este simple analisis, podemos ver que la temperatura media y su cuadrado tiene el mayor efecto, seguido de la lluvia y las rachas de viento. Quiero ver como de preciso son capaces de adivinar estas variables la demanda real. Lo normal es que no obtengamos un buen resultado, ya que solo contenemos 357 dias de observaciones. Por lo general, con series temporales se suelen realizar estos analisis con 5 anos de datos.

#Realizamos un grafico Demanda Real vs. Demanda predicha
ggplot(datos_modelo,aes(x=Fitted, y=Demanda))+geom_point()+xlab("Demanda predicha")+ylab("Demanda Real")+geom_abline()+ggtitle("Observaciones Demanda Real vs. Predicha")

#Realizamos un grafico de los valores predecidos
datos_modelo$Date<-as.Date(datos_modelo$Date)
ggplot(datos_modelo, aes(x=Date)) + 
  geom_line(aes(y = Demanda), color = "red") + 
  geom_line(aes(y = Fitted), color="blue", linetype="twodash") +ggtitle("Valores Demanda Real(rojo) vs Valores Demanda Predicha (azul)")

#Modelos no lineales - en vista al analisis descriptivo del clima, podriamos intentar construir un modelo no lineal con una funcion S-curve, aunque no mejoraria demasiado. Ademas podemos concluir que este modelo no serviria para machine learning (muy impreciso) debido a los altos residuos en esta muestra y un coeficiente de R^2 que apenas explica un 41% de la variable dependiente.

#scurve = function(x, center, width) {1 / (1 + exp(-(x - center) / width))}

11-Analisis del clima y demanda diaria

Analisis de la variable Demanda Diaria respecto a diferentes tipos de variables climaticas

#Temperatura media respecto a demanda diaria
ggplot(datos_modelo, aes(Tmean,Demanda)) + geom_smooth() + labs(x = "Temperatura Media (centigrados)", y = "Demanda diaria")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#Racha de viento respecto a demanda diaria
ggplot(datos_modelo, aes(Wind_Gust,Demanda)) + geom_smooth() + labs(x = "Racha del viento (m/s)", y = "Demanda diaria")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#Velocidad del viento respecto a demanda diaria
ggplot(datos_modelo, aes(Wind_speed,Demanda)) + geom_smooth() + labs(x = "Velocidad del viento (m/s)", y = "Demanda diaria")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 17 rows containing non-finite values (stat_smooth).

#Precipitaciones respecto a demanda diaria
ggplot(datos_modelo, aes(Rain_amount,Demanda)) + geom_smooth() + labs(x = "Precipitacion por metro cuadrado (mm = l/m2)", y = "Demanda diaria")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#Temperatura media respecto al mes
datos_modelo$Month<-as.factor(datos_modelo$Month)
ggplot(datos_modelo, aes(Month,Tmean)) + geom_violin() + labs(y = "Temperatura Media (centigrados)", x = "Mes")

#Temperatura maxima respecto a Demanda diaria
ggplot(datos_modelo, aes(Tmax,Demanda)) + geom_smooth() + labs(y = "Demanda diaria", x = "Temperatura Maxima (centigrados)")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

rm(datos_clima)

12-Movimientos por hora

En esta seccion intentaremos explicar a que horas se mueven mas los usuarios y a que hora se mueven mas los empleados

NBicis_usuarios<-count(BiciMAD_paseyocasional,unplug_hourTime)
NBicis_empleados<-count(BiciMAD_empleados,unplug_hourTime)

NBicis_usuarios<-setNames(NBicis_usuarios,c("Hora","Numero_usuarios"))
NBicis_empleados<-setNames(NBicis_empleados,c("Hora","Numero_empleados"))

Uso_hora<-merge(NBicis_empleados,NBicis_usuarios, by="Hora")

Uso_hora$Hora<-as.factor(Uso_hora$Hora)
Uso_hora$Numero_empleados <-as.numeric(as.character(Uso_hora$Numero_empleados))
Uso_hora$Numero_usuarios<-as.numeric(as.character(Uso_hora$Numero_usuarios))

ggplot(Uso_hora)+geom_line(aes(Hora,Numero_usuarios, group=1, col="Usuarios"))+geom_line(aes(Hora,Numero_empleados, col="Empleados", group=2))+ylab("Numero de Bicicletas")+xlab("Hora")+scale_y_continuous(labels = comma)+ggtitle("Numero de usuarios (azul) vs. Numero de empleados (rojo) por hora")

13-Modelos no supervisados

Realizamos un analisis por clusters despues de multiples transformaciones de los datos facilitados por el ayuntamiento de Madrid. Se obtuvo un registro en formato JSON, que se transformo y formateo a XLSX, con las bicicletas disponibles, bases libres donde anclar, numero de reservas y mas, por dia y por hora del mes de septiembre. Se computaron las medias a traves de una tabla dinamica. Estas medias son las usadas para analizar las estaciones por clusters.

K-means Clustering

13.1-Primera agrupacion por clusters

Agrupacion por medias mensuales, plano de 8 dimensiones (4 semanas, 2 variables)

estaciones_sep<-read.csv("estaciones_septiembre_mensual.csv", sep=";")
str(estaciones_sep)
## 'data.frame':    172 obs. of  10 variables:
##  $ nombre.estacion       : Factor w/ 172 levels "Agustin de Betancourt",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ avg.ocupacion.ancladas: num  0.298 0.424 0.462 0.234 0.357 ...
##  $ avg.bases.libres      : num  0.702 0.576 0.536 0.766 0.643 ...
##  $ avg.capacidad.estacion: num  0.981 0.984 0.983 0.97 0.973 ...
##  $ avg.bases.inoperativas: num  0.449 0.381 0.452 0.713 0.659 ...
##  $ avg.bases.operativas  : num  23.6 23.6 26.5 23.3 23.3 ...
##  $ avg.bases.totales     : num  23.9 23.9 24.5 23.8 23.8 ...
##  $ longitud.estacion     : num  -3.7 -3.68 -3.68 -3.67 -3.69 ...
##  $ lat.estacion          : num  40.4 40.5 40.4 40.4 40.4 ...
##  $ avg.numero.reservas   : num  0.1057 0.0751 0.057 0.0223 0.0974 ...
#Cargamos la libreria
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
#Pondremos la variable del nombre de la estacion como nombre de columna
estaciones_sep_rows <- data.frame(estaciones_sep[,-1], row.names=estaciones_sep[,1])

#Seleccionamos las variables
estaciones_sep_rows$lat.estacion<-NULL
estaciones_sep_rows$longitud.estacion<-NULL
#estaciones_sep_rows_var <- estaciones_sep_rows[,1:28]

#Quitaremos las estaciones no operativas, Red de San Luis, Carretas y Anton Martin
estaciones_sep_rows<-estaciones_sep_rows[!(estaciones_sep_rows$avg.ocupacion.ancladas==0.0000000 & estaciones_sep_rows$avg.bases.libres==0.0000000),]


#Escalamos las variables seleccionadas del dataframe a las mismas dimensiones
estaciones_sep_scaled <- scale(estaciones_sep_rows)


#Obtenemos las distancias:
distancia <- get_dist(estaciones_sep_scaled)
fviz_dist(distancia,lab_size = 5.2)

# Encontramos el numero de clusters: K optimo
fviz_nbclust(estaciones_sep_scaled,kmeans,method="wss") # Elbow

fviz_nbclust(estaciones_sep_scaled,kmeans,method="silhouette") # Silhouette

# Kmeans
a <- kmeans(estaciones_sep_scaled , 5 , nstart=15)
a$cluster
##                         Agustin de Betancourt 
##                                             5 
##                               Alberto Alcocer 
##                                             5 
##                                        Alcala 
##                                             2 
##                                     Alcantara 
##                                             5 
##                                       Almaden 
##                                             5 
##                               Alonso Martinez 
##                                             5 
##                                    Altamirano 
##                                             1 
##                                 Antonio Maura 
##                                             5 
##                          Arcipreste de Hita A 
##                                             5 
##                          Arcipreste de Hita B 
##                                             5 
##                                      Atocha A 
##                                             1 
##                                      Atocha B 
##                                             2 
##                              Augusto Figueroa 
##                                             5 
##                                         Ayala 
##                                             5 
##                             Banco de Espana A 
##                                             2 
##                                       Barcelo 
##                                             5 
##                           Biblioteca Nacional 
##                                             5 
##                              Bravo Murillo 44 
##                                             5 
##                                   Calle Mayor 
##                                             2 
##                             Carlos Cambronero 
##                                             3 
##                                Casa Encendida 
##                                             1 
##                                    Castellana 
##                                             5 
##                                Castellana 106 
##                                             5 
##                                Castellana 164 
##                                             1 
##                                 Castellana 42 
##                                             3 
##           Castellana frente a Hermanos Pinzon 
##                                             5 
##                                      Castello 
##                                             2 
##                                  Cea Bermudez 
##                                             2 
##                                       Cibeles 
##                                             3 
##                                Claudio Coello 
##                                             5 
##                           Colegio Arquitectos 
##                                             5 
##                                       Colon A 
##                                             4 
##                                       Colon B 
##                                             4 
##                                      Colombia 
##                                             1 
##                                      Columela 
##                                             5 
##                                Conde de Casal 
##                                             1 
##                            Conde de Romanones 
##                                             1 
##                                   Conde Duque 
##                                             5 
##                                Conde Penalver 
##                                             5 
##                              Cuesta de Moyano 
##                                             5 
##                              Diego de Leon 52 
##                                             3 
##                               Doce de Octubre 
##                                             1 
##                                Doctor Arce 45 
##                                             5 
##                            Embajada de Italia 
##                                             5 
##                                 Embajadores 1 
##                                             1 
##                                 Embajadores 2 
##                                             1 
##                              Entrada Matadero 
##                                             1 
##                                      Espalter 
##                                             5 
##                          Fernando el Catolico 
##                                             1 
##                                        Ferraz 
##                                             5 
##                                    Fuencarral 
##                                             2 
##                     General Alvarez de Castro 
##                                             1 
##                      General Martinez Campos  
##                                             5 
##                              General Pardinas 
##                                             2 
##                               General Peron 1 
##                                             1 
##         General Peron con Poeta Joan Maragall 
##                                             5 
##                Glorieta de los Cuatro Caminos 
##                                             1 
##                               Guzman el Bueno 
##                                             5 
##                                     Hortaleza 
##                                             3 
##                                       Huertas 
##                                             3 
##                                         Ibiza 
##                                             1 
##                             Jacinto Benavente 
##                                             5 
##                                   Jacometrezo 
##                                             3 
##                                 Jesus y Maria 
##                                             1 
##                                    Jorge Juan 
##                                             1 
##                                  Jose Abascal 
##                                             1 
##                       Jose Gutierrez Abascal  
##                                             5 
##                                 Juan Bravo 50 
##                                             1 
##                                   Juan Martin 
##                                             1 
##                                      Malasana 
##                                             3 
##                                Manuel Silvela 
##                                             1 
##                             Maria Francisca 1 
##                                             4 
##                              Marques de Cubas 
##                                             3 
##                        Marques de la Ensenada 
##                                             5 
##                          Marques de Salamanca 
##                                             5 
##                               Menendez Pelayo 
##                                             1 
##                               Metro Piramides 
##                                             1 
##                                 Mendez Alvaro 
##                                             1 
##                      Martires Concepcionistas 
##                                             2 
##                                   Miguel Moya 
##                                             5 
##                                       Moncloa 
##                                             5 
##                             Museo Reina Sofia 
##                                             5 
##                                       Narvaez 
##                                             1 
##                                    Nino Jesus 
##                                             1 
##                                     O'Donnell 
##                                             1 
##                                     Orense 12 
##                                             5 
##                                     Orense 36 
##                                             5 
##                               Ortega y Gasset 
##                                             3 
##                            Ortega y Gasset 87 
##                                             1 
##                           Palacio de Deportes 
##                                             5 
##                            Palacio de Oriente 
##                                             5 
##                          Palos de la Frontera 
##                                             1 
## Paseo de la Castellana con Raimundo Fernandez 
##                                             5 
##                         Paseo de la Esperanza 
##                                             1 
##                           Paseo de la Florida 
##                                             1 
##                         Paseo de la Habana 42 
##                                             5 
##                         Paseo de la Habana 63 
##                                             5 
##                         Paseo de las Delicias 
##                                             1 
##                                Paseo de Moret 
##                                             5 
##                                   Pedro Bosch 
##                                             1 
##      Castellana - Glorieta de Emilio Castelar 
##                                             4 
##                                    Pio Baroja 
##                                             5 
##                                Pintor Rosales 
##                                             3 
##                            Plaza Conde Suchil 
##                                             4 
##                           Plaza de Celenque A 
##                                             5 
##                           Plaza de Celenque B 
##                                             5 
##                                Plaza de Colon 
##                                             5 
##                             Plaza de Espana A 
##                                             5 
##                             Plaza de Espana B 
##                                             5 
##                            Plaza de Felipe II 
##                                             5 
##                           Plaza de Juan Pujol 
##                                             3 
##                            Plaza de la Cebada 
##                                             1 
##                     Plaza de la Independencia 
##                                             5 
##                         Plaza de la Provincia 
##                                             4 
##               Plaza de la Republica Argentina 
##                                             5 
##                          Plaza de las Salesas 
##                                             3 
##                             Plaza de Lavapies 
##                                             1 
##                           Plaza de los Astros 
##                                             1 
##                           Plaza de los Carros 
##                                             1 
##                       Plaza de Nelson Mandela 
##                                             1 
##                         Plaza de Pedro Zerolo 
##                                             3 
##                              Plaza de Ramales 
##                                             5 
##                        Plaza de San Ildefonso 
##                                             3 
##                  Plaza de San Juan de la Cruz 
##                                             5 
##                           Plaza de San Miguel 
##                                             5 
##                            Plaza de Santa Ana 
##                                             3 
##                             Plaza del Carmen  
##                                             5 
##                              Plaza del Cordon 
##                                             5 
##                           Plaza San Francisco 
##                                             1 
##                           Principe de Vergara 
##                                             5 
##                                          Prim 
##                                             3 
##                             Puerta de Granada 
##                                             1 
##                              Puerta de Madrid 
##                                             2 
##                    Puerta de Mariano de Cavia 
##                                             1 
##                              Puerta de Toledo 
##                                             1 
##                    Puerta del Doce de Octubre 
##                                             1 
##                        Puerta del Angel Caido 
##                                             2 
##                          Puerta del Pacifico  
##                                             1 
##                              Puerta del Sol A 
##                                             5 
##                              Puerta del Sol B 
##                                             1 
##                                      Quintana 
##                                             1 
##                            Raimundo Fernandez 
##                                             5 
##                          Ribera de Curtidores 
##                                             1 
##                               Ronda de Atocha 
##                                             2 
##                              Sainz de Baranda 
##                                             1 
##                                    San Andres 
##                                             5 
##                                  San Bernardo 
##                                             5 
##                                    San German 
##                                             1 
##                              San Hermenegildo 
##                                             3 
##                                   San Quintin 
##                                             2 
##                            San Vicente Ferrer 
##                                             3 
##                      Santa Cruz del Marcenado 
##                                             5 
##                            Santa Engracia 127 
##                                             4 
##                             Santa Engracia 14 
##                                             5 
##                                  Santa Isabel 
##                                             1 
##                      Santa Maria de la Cabeza 
##                                             1 
##                                 Santo Domingo 
##                                             3 
##                                    Segovia 26 
##                                             5 
##                                    Segovia 45 
##                                             1 
##                                       Serrano 
##                                             5 
##                                   Serrano 210 
##                                             5 
##                                       Sevilla 
##                                             3 
##                         Sor Angela de la Cruz 
##                                             5 
##                                   Valderribas 
##                                             1 
##                                     Velazquez 
##                                             5 
##                                 Velazquez 130 
##                                             1 
##                             Ventura Rodriguez 
##                                             5 
##                                    Villanueva 
##                                             5
a$size
## [1] 55 13 21  7 72
a$centers
##   avg.ocupacion.ancladas avg.bases.libres avg.capacidad.estacion
## 1             1.11416960      -1.10823906             0.24715213
## 2            -0.08132129       0.09117698            -0.05980603
## 3            -0.86070261       0.80462881            -1.99949743
## 4             0.24085109      -0.24273862            -0.21215870
## 5            -0.60879659       0.61902518             0.42581484
##   avg.bases.inoperativas avg.bases.operativas avg.bases.totales
## 1             -0.2379264            0.1369591        0.06029540
## 2              0.3102350            1.6951475        1.83411536
## 3              1.9649855           -0.8117377       -0.20312170
## 4             -0.2384054           -3.1206066       -3.24523850
## 5             -0.4242078            0.1294593       -0.00246558
##   avg.numero.reservas
## 1          0.68113291
## 2         -0.05519189
## 3         -0.02662197
## 4         -0.34738428
## 5         -0.46880644
fviz_cluster(a, data=estaciones_sep_scaled, labelsize = 7)

# Agrego cluster asociado
estaciones_sep_rows$cluster <- a$cluster

# Resumen por cluster

estaciones_sep<-estaciones_sep[!(estaciones_sep$avg.ocupacion.ancladas==0.0000000 & estaciones_sep$avg.bases.libres==0.0000000),]

aggregate(estaciones_sep,list(estaciones_sep_rows$cluster),mean)
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
##   Group.1 nombre.estacion avg.ocupacion.ancladas avg.bases.libres
## 1       1              NA              0.5338822        0.4657891
## 2       2              NA              0.4053700        0.5944160
## 3       3              NA              0.3215886        0.6709275
## 4       4              NA              0.4400027        0.5586065
## 5       5              NA              0.3486678        0.6510231
##   avg.capacidad.estacion avg.bases.inoperativas avg.bases.operativas
## 1              0.9651386              0.8334050             23.09124
## 2              0.9569099              1.1806997             26.19835
## 3              0.9049125              2.2290880             21.19948
## 4              0.9528258              0.8331015             16.59547
## 5              0.9699280              0.7153840             23.07628
##   avg.bases.totales longitud.estacion lat.estacion avg.numero.reservas
## 1          23.82066         -3.631824     39.75582          0.12239221
## 2          24.63464         -3.692348     40.42044          0.08484006
## 3          23.69978         -3.540992     38.69108          0.08629711
## 4          22.30379         -3.694011     40.43093          0.06993841
## 5          23.79186         -3.696223     40.42950          0.06374594
a$size
## [1] 55 13 21  7 72

5 grupos de estaciones con comportamientos operativos muy diferentes entre sí.

Hierarchical Clustering

13.2-Segunda agrupacion por clusters

Utilizaremos los datos mensuales del punto 13.1. Pero esta vez utilizamos clustering jerarquico

#Calculamos las distancias
static.hc <- hclust(dist(estaciones_sep_scaled))
fviz_dist(dist(estaciones_sep_scaled), lab_size = 4)

#Como hemos visto hay algunas estaciones que han saltado del grupo 2 al 1. Eso significa que existen una especie de estaciones entre medias. Decidimos construir el dendograma con 3 clusters.
static.hc <- hcut(estaciones_sep_scaled, k = 5)

#Pintamos el dendrograma
fviz_dend(static.hc, rect = TRUE, cex = 0.5, k_colors = c("#E69F00", "#56B4E9","blue","black", "red"))

14-Mapa de Clusters

Vamos a realizar ahora un mapa con las estaciones por cluster, para ver si podemos ver algun patron a nivel geografico.

#Sacamos la categoria de cluster de nuestra variable a
cluster_list<-as.data.frame(a$cluster)
cluster_list<-setNames(cluster_list,"cluster")

cluster_listhc <-as.data.frame(static.hc$cluster)
cluster_listhc<-setNames(cluster_listhc,"clusterhc")

#Convertimos las row.names a una columna con el nombre de la estacion
cluster_list<-as_tibble(rownames_to_column(cluster_list, var = "nombre.estacion"))

cluster_listhc<-as_tibble(rownames_to_column(cluster_listhc, var = "nombre.estacion"))

#Lo unimos a nuestro dataset original
#estaciones_sep<-setNames(estaciones_sep, c("nombre.estacion","avg.ocupacion.ancladas","avg.bases.libres","avg.capacidad.estacion", "avg.bases.inoperativas","avg.bases.operativas" etc.))


cluster_data<-left_join(cluster_list,cluster_listhc, by="nombre.estacion")

rm(cluster_listhc)
rm(cluster_list)


cluster_dataset<-left_join(cluster_data,estaciones_sep, by="nombre.estacion")

#Crearemos un mapa con la nueva columna

#Sacamos el mapa de Madrid de Google Maps con nuestra API y ggmap
mapa_clusters<-ggmap(get_googlemap(center = c(lon = -3.7, lat = 40.43), zoom = 13, scale = 2, maptype ="terrain", color = "bw")) 
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.43,-3.7&zoom=13&size=640x640&scale=2&maptype=terrain&key=xxx
#Crearemos un mapa con los clusters
cluster_dataset$cluster<-as.factor(cluster_dataset$cluster)
cluster_dataset$clusterhc<-as.factor(cluster_dataset$clusterhc)
cluster_dataset$longitud.estacion<-as.numeric(as.character(cluster_dataset$longitud.estacion))
cluster_dataset$lat.estacion<-as.numeric(as.character(cluster_dataset$lat.estacion))

mapa_clusters + geom_point(aes(y =lat.estacion , x = longitud.estacion , color = cluster), data=cluster_dataset)+scale_color_manual(values=c("#E69F00", "#56B4E9","blue","black", "red"))
## Warning: Removed 2 rows containing missing values (geom_point).

mapa_clusters + geom_point(aes(y =lat.estacion , x = longitud.estacion , color = clusterhc), data=cluster_dataset)+scale_color_manual(values=c("#E69F00", "#56B4E9","blue","black", "red"))
## Warning: Removed 2 rows containing missing values (geom_point).

Como vemos las zonas mas alejadas y algunas proximas al centro son las de un cluster. Toda la zona sur que esta muy proxima entre si puede servir para reponer las biciletas de la zona centro, donde esta la mayoria de los demas clusteres. El mayor problema esta al norte del Parque de El Retiro, donde encontramos poca proximidad entre estaciones con alto numero de ocupacion, especialmente en el barrio de Salamanca. No obstante, al este del Retiro existe una zona donde todas las estaciones gozan de disponibilidad y siempre hay bicicletas. Ademas, se podria aprovechar la vertical de Dr.Esquerdo para subir estas bicicletas a la zona de Salamanca.

Nota final

Cabe destacar que solo se ha analizado el periodo de 2018, ya que 2017 y 2019 se muestran incompletos en el portal OPENDATA. Al tratar con series temporales el comportamiento de los usuarios no es el mismo (e.j en Invierno hay menos viajes y posiblemente a estaciones mas cercanas debido al frio). Por ello, al disponer unicamente de este ano al completo, se ha elaborado el estudio con estos mismos datos. En 2019, se ha realizado la primera expansion de BiciMAD hacia Tetuan como exponemos en las observaciones, por lo que habria que realizar este mismo estudio cuando salgan los datos. En 2020, se ha acordado la siguiente expansion del sistema mas alla de la M-30.

El analisis continua en el rebalanceo por clustering con series temporales.

Noticias recientes sobre BiciMAD: https://elpais.com/ccaa/2019/10/13/madrid/1570998467_632903.html